#####################################
# Replication for municipality-level analysis for 
# Charnysh and Schaub, Migration and Social Change: Evidence from post-WWII Displacement in Germany, JOP
# R version 4.3.1 (2023-06-16) 
#####################################

rm(list = ls())

# Load packages
load <- c("readstata13", "openxlsx", "sandwich", "lmtest", "stargazer", "dplyr", "conleyreg","sensemakr", "reshape2", "DirectEffects", "texreg","xtable", "interflex", "haven", "sf", "ggplot2", "geosphere", "reldist", "ape", "GWmodel")
lapply(load, library, character.only = TRUE)


####### CREATE FUNCTIONS ########
cse = function(reg) {
  rob = sqrt(diag(vcovHC(reg, type = "HC1")))
  return(rob)
}


#function to create formulas
formulas <- function(a, b) {
  formulas <- list()
  n <- 0
    for (i in a) {
      for (j in b) {
        n <- n + 1
        formulas[n] <- paste(i, j)
      }
    }
  return(formulas)
}


#function to run models
run_models <- function(formlist) {
  models <- list()
  n <- 0
  for (i in 1:length(formlist)){
    n <- n + 1
    models[[n]] <- eval(parse(text=formlist[[n]])) 
  }
  return(models)
}

#vcov function 
hc1<-function(x) vcovHC(x, type="HC1")  

#function to extract coefficients + calculate errors
run_coef <- function(output) {
  stdev <- list() 
  n <- 0
  for (i in 1:length(output)){
    n <- n + 1
    stdev[[n]] <- coeftest(output[[n]], vcov. = hc1)
  }
  return(stdev)    
}

#Function to format pvalues as stars for output from sequential g-estimation
pval_stars <- function(p) {
  if (p < 0.01) {
    return("***")
  } else if (p < 0.05) {
    return("**")
  } else if (p < 0.1) {
    return("*")
  } else {
    return("")
  }
}
########################################
######## Compute main variables: #######
########################################

dat<-read.xlsx("municipal_data/municipality-dataset-short.xlsx")

dat$sprotestants1939<-dat$protestants1939/dat$pop1939orig #share of protestants in 1939
dat$scatholics1939<-dat$catholics1939/dat$pop1939orig #share of catholics in 1939
dat$sotherrel1939<-(dat$pop1939orig - (dat$catholics1939 + dat$protestants1939))/dat$pop1939orig #share of other religious groups in 1939
dat$relfrac1939<-1-(dat$sprotestants1939^2+dat$scatholics1939^2+dat$sotherrel1939^2) #religious fractionalization in 1939

dat$sprotestants1950<-dat$protestants1950/dat$pop1950 
dat$scatholics1950<-dat$catholics1950/dat$pop1950
dat$sotherrel1950<-(dat$pop1950 - (dat$catholics1950 + dat$protestants1950))/dat$pop1950
dat$relfrac1950<-1-(dat$sprotestants1950^2+dat$scatholics1950^2+dat$sotherrel1950^2)

dat$DeltaFrac<-dat$relfrac1950-dat$relfrac1939

#Note we use pop1939orig from 1939 census publication for religion to ensure that shares add up to 1.

#### Main paper:
dat$lnpop1939<-log(dat$pop1939) 
dat$lnpop1950<-log(dat$pop1950)

#intermediate confounders for g-estimation
dat$PT65andolder1950<-dat$aged65up1950/dat$pop1950

dat$ShareCommuters1950<-(dat$commuters_ex1950+dat$commuters_in1950)/dat$pop1950 #in & out commuters
dat$taxRevenuePC1950<-dat$taxrevenues1950/dat$pop1950

dat$OtherReligion_1970<-dat$Bevoelkerung_Religion_insg_1970-dat$Katholisch_1970-dat$Evangelisch_1970

dat$relfrac1970<-1-((dat$Katholisch_1970/dat$Bevoelkerung_Religion_insg_1970)^2+(dat$Evangelisch_1970/dat$Bevoelkerung_Religion_insg_1970)^2+(dat$OtherReligion_1970/dat$Bevoelkerung_Religion_insg_1970)^2)
dat$relfrac1970<-ifelse(dat$relfrac1970<0, NA, dat$relfrac1970) #Replace number that is too high and must be a typo with NA. 

dat$lndistBorderEast <-log(dat$distBorderEast)

dat$srefugees1950<-dat$prefugees1950/100 #express percentage as share for consistency with other variables

dat$InAgric_50<-dat$wagriforestry1950/dat$nactive1950 
dat$InAgric_70<-dat$LForstwirtschaft_insg_1970/dat$Erwerbstaetige_insg_1970

dat$InManuf70 <-dat$Produzierenes_Gewerbe_insg_1970/dat$Erwerbstaetige_insg_1970

dat$HandelVerkehr70 <-dat$Handel_Verkehr_insg_1970/dat$Erwerbstaetige_insg_1970

dat$pAuslaender1970<-dat$Auslaender_insg_1970/dat$Bevoelkerung_insges_1970

dat$SelbstStandige_70<-dat$Selbstaendige_insg_1970/dat$Erwerbstaetige_insg_1970


dat$attendants_rel1970s<-dat$attendants_rel1970/100 #express percentage as share
dat$attendants_rel1980s<-dat$attendants_rel1980/100
dat$attendants_rel1990s<-dat$attendants_rel1990/100

dat$preligious_50s<-dat$preligious_50/100 #express percentage as share
dat$preligious_70s<-dat$preligious_70/100
dat$preligious_87s<-dat$preligious_87/100

dat$pdivorce_tomar_70<-dat$Geschieden_insg_1970/dat$Verheiratet_insg_1970
dat$pdivorce_tomar_87<-dat$Geschieden_insg_1987/dat$Verheiratet_insg_1987

#calculate party vote shares to eligible voters for CSU:
dat$csu_voters1949<-dat$votes_csu_1949/dat$eligible_voters_1949
dat$csu_voters1953<-dat$votes_csu_1953/dat$eligible_voters_1953  
dat$csu_voters1957<-dat$votes_csu_1957/dat$eligible_voters_1957
dat$csu_voters1961<-dat$votes_csu_1961/dat$eligible_voters_1961 
dat$csu_voters1965<-dat$votes_csu_1965/dat$eligible_voters_1965  
dat$csu_voters1969<-dat$votes_csu_1969/dat$eligible_voters_1969 
dat$csu_voters1972<-dat$votes_csu_1972/dat$eligible_voters_1972  
dat$csu_voters1976<-dat$votes_csu_1976/dat$eligible_voters_1976 
dat$csu_voters1980<-dat$votes_csu_1980/dat$eligible_voters_1980 
dat$csu_voters1983<-dat$votes_csu_1983/dat$eligible_voters_1983
dat$csu_voters1987<-dat$votes_csu_1987/dat$eligible_voters_1987 
dat$csu_voters1990<-dat$votes_csu_1990/dat$eligible_voters_1990

#calculate party vote shares to eligible voters for the Green Party:
dat$greens_voters1980<-dat$votes_greens_1980/dat$eligible_voters_1980 
dat$greens_voters1983<-dat$votes_greens_1983/dat$eligible_voters_1983
dat$greens_voters1987<-dat$votes_greens_1987/dat$eligible_voters_1987 
dat$greens_voters1990<-dat$votes_greens_1990/dat$eligible_voters_1990


#Create factor for fixed effects used in Conley error specifications 
dat$county_id_f<-as.factor(dat$kreis_kenn1950_alt) #prewar kreise had same borders as 1950 kreise
dat$diocese_id_f<-as.factor(dat$diocese) #diocese fixed effect


###################################
######## FIGURE 2 in the paper ####
###################################

#jpg("Change_fractionalization.jpg", width = 800, height = 600, units = "px", res = 300)
tiff("Figure2.tiff", width = 8, height = 6, units = "in", res = 300)
plot(dat$relfrac1939, dat$relfrac1950, 
     ylab="Religious fractionalization in 1950", 
     xlab="Religious fractionalization in 1939", 
     pch=20)
abline(0, 1, col="red", lwd=3)
dev.off()

######################################################
###########  TABLE A1: Descriptive statistics  #######
######################################################

vars<-c("DeltaFrac", "srefugees1950", "relfrac1950",  "relfrac1939",  "scatholics1950",  "scatholics1939","sprotestants1950",  "sprotestants1939", "sfemale1939",  "sagri1939", "diststation1950",  "pop1939", "pop1950",   "distBorderEast", "dest1945", "popdens1939", "InAgric_50")

dvs= c("attendants_rel1970s", "attendants_rel1980s", "attendants_rel1990s", "preligious_70s", "preligious_87s", "pdivorce_tomar_70", "pdivorce_tomar_87")


stargazer(dat[, vars], title="Descriptive statistics for Bavarian municipalities", digits=3, summary.stat=c("n", "mean","sd",  "min", "max"),  label = "tab:descriptive",
          covariate.labels=c("Delta Diversity (1939 to 1950)", "Share expellees (1950)","Religious Fractionalization (1950)", "Religious Fractionalization (1939)", 
                             "Share Catholics (1950)", "Share Catholics (1939)","Share Protestants (1950)", "Share Protestants (1939)", "Share female (1939)", 
                             "Share in agriculture (1939)", "Distance to station (1950)",  "Population (1939)", "Population (1950)", "Distance to Eastern Border", "Destruction (1945)",
                             "Population Density (1939)","Share in agriculture (1950)"))

#This table is not in the paper because we show DV mean and sd in regression tables:
stargazer(dat[, dvs], title="Descriptive statistics for Bavarian municipalities", digits=3, summary.stat=c("n", "mean","sd",  "min", "max"),  label = "tab:descriptive_dv",
          covariate.labels=c("Mass attendance (1970)", "Mass attendance (1980)","Mass attendance (1990)", "Church affiliation (1970)", 
                             "Church affiliation (1980)", "Church affiliation (1987)", "Divorces to marriages (1970)", "Divorces to marriages (1987)"))





#########################################################################
###### Table A3: Correlations and variance explained ####################
#########################################################################


munich_coords <- c(11.5820, 48.1351)  # Munich (Longitude, Latitude)
dat$dist_to_munich <- distHaversine(
  matrix(c(dat$lon, dat$lat), ncol = 2),
  munich_coords
) / 1000


dat$lndist_to_munich<-log(dat$dist_to_munich)
dat$serbhofe1939<-dat$nerbhofe1939/dat$nfarms1939
dat$PtFarmsUnder2ha<-dat$nfarms05_2ha1939/dat$nfarms1939
dat$PtFarms2to5ha<-dat$nfarms2bis5ha1939/dat$nfarms1939
dat$PtFarms5to10ha<-dat$nfarms5bis10ha1939/dat$nfarms1939
dat$PtFarms10to20ha<-dat$nfarms10bis20ha1939/dat$nfarms1939
dat$PtFarms20to100ha<-dat$nfarms20bis100ha1939/dat$nfarms1939
dat$PtFarmsOver100ha<-dat$nfarms100ha1939/dat$nfarms1939

dat$UrbanType<-ifelse(dat$Type=="Landkreis", 0, 1)
dat$land.gini1939 <- apply(dat, 1, function(row) {
  gini(x = c(1, 3.5, 7.5, 15, 60, 120), 
       weights = as.numeric(row[c("PtFarmsUnder2ha", "PtFarms2to5ha", "PtFarms5to10ha",
                                  "PtFarms10to20ha", "PtFarms20to100ha", "PtFarmsOver100ha")]))
})

# subset to non-missing
dat_res<-subset(dat, !is.na(relfrac1939) & !is.na(serbhofe1939))


model <- lm(relfrac1950 ~ relfrac1939 +lndistBorderEast+dest1945, data = dat_res)
residuals <- residuals(model)

dat_res$residuals<-residuals(model)

corrsResiduals<-rbind(cbind(cor.test(dat_res$residuals, dat_res$sfemale1939)$estimate, summary(lm(residuals~ sfemale1939, data=dat_res))$r.squared), 
                      cbind(cor.test(dat_res$residuals, dat_res$sagri1939)$estimate,summary(lm(residuals~ sagri1939, data=dat_res))$r.squared), 
                      cbind(cor.test(dat_res$residuals, dat_res$land.gini1939)$estimate, summary(lm(residuals~ land.gini1939, data=dat_res))$r.squared),
                      cbind(cor.test(dat_res$residuals, dat_res$serbhofe1939)$estimate, summary(lm(residuals~ serbhofe1939, data=dat_res))$r.squared),
                      cbind(cor.test(dat_res$residuals, dat_res$lndist_to_munich)$estimate, summary(lm(residuals~ lndist_to_munich, data=dat_res))$r.squared),
                      cbind(cor.test(dat_res$residuals, dat_res$diststation1950)$estimate, summary(lm(residuals~ diststation1950, data=dat_res))$r.squared),
                      cbind(cor.test(dat_res$residuals, dat_res$UrbanType)$estimate, summary(lm(residuals~ UrbanType, data=dat_res))$r.squared),
                      cbind(cor.test(dat_res$residuals, dat_res$pop1939)$estimate, summary(lm(residuals~ pop1939, data=dat_res))$r.squared), #not sig
                      cbind(cor.test(dat_res$residuals, dat_res$lnpop1939)$estimate, summary(lm(residuals~ lnpop1939, data=dat_res))$r.squared), #not sig.
                      cbind(cor.test(dat_res$residuals, dat_res$popdens1939)$estimate, summary(lm(residuals~ popdens1939, data=dat_res))$r.squared),
                      cbind(cor.test(dat_res$residuals, dat_res$altitude1950)$estimate, summary(lm(residuals~ altitude1950, data=dat_res))$r.squared),
                      cbind(cor.test(dat_res$residuals, dat_res$lat)$estimate, summary(lm(residuals~ lat, data=dat_res))$r.squared),
                      cbind(cor.test(dat_res$residuals, dat_res$lon)$estimate, summary(lm(residuals~ lon, data=dat_res))$r.squared))

corrsResiduals<-as.data.frame(corrsResiduals, 
                              row.names=c("Share female (1939)","Share in agriculture (1939)",  "Landholding gini (1939)", "Share of hereditary farms (1939)",
                                           "Ln distance to Munich", "Distance to station (1950)", "City (Stadtkreis)", "Population (1939)", "Ln population (1939)", "Pop density (1939)", "Altitude","Latitude", "Longitude")) 

names(corrsResiduals)<-c("Pearson Correlation", "Propotion of explained variance (R^2)")

xtable(corrsResiduals, digits=3)

#########################################################################
###### Figure A2 (top): Moran’s I statistic for the full model ##########
#########################################################################

coords <- cbind(dat_res$lon, dat_res$lat)
distance_matrix <- gw.dist(dp.locat = coords, longlat = TRUE)
threshold <- 15  # 15km threshold 
weight_matrix <- gw.dist(dp.locat = coords, longlat = TRUE)
weight_matrix[weight_matrix > threshold] <- 0 
weight_matrix[weight_matrix > 0] <- 1  

# Residualization using controls from the paper.
model <- lm(relfrac1950 ~ relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data = dat_res)
dat_res$residuals_full<-residuals(model)

## Analytic results for Moran's I, non-residualized vs. residualized values
Moran.I(dat_res$DeltaFrac, weight_matrix)
Moran.I(dat_res$residuals_full, weight_matrix)


## Randomization inference

# Seed and N simulations
set.seed(7) 
n_simulations <- 1000


# Function for randomization inference for Moran's I. Note it takes some time to run.
run_moran_randomization <- function(data, exp_weighted_matrix, n_simulations) {
  randomized_moran <- numeric(n_simulations)
  for (i in seq_len(n_simulations)) {
    shuffled_data <- sample(data)
    moran_result <- Moran.I(shuffled_data, exp_weighted_matrix)
    randomized_moran[i] <- moran_result$observed
  }
  return(randomized_moran)
}

# Values needed for plotting
moran_residualized <- Moran.I(dat_res$residuals_full, weight_matrix)
obs_stat_moran_res <- moran_residualized$observed
dat_res$DeltaFrac<-dat_res$relfrac1950-dat_res$relfrac1939
moran_non_residualized <- Moran.I(dat_res$DeltaFrac, weight_matrix)
obs_stat_moran_nr <- moran_non_residualized$observed

randomized_moran_residualized <- run_moran_randomization(dat_res$residuals_full, weight_matrix, n_simulations)
randomized_moran_non_residualized <- run_moran_randomization(dat_res$DeltaFrac, weight_matrix, n_simulations)

# Data frame for plotting
df_moran_residualized <- data.frame(stat = randomized_moran_residualized, group = "Residualized")
df_moran_non_residualized <- data.frame(stat = randomized_moran_non_residualized, group = "Non-Residualized")
df_moran_combined <- rbind(df_moran_residualized, df_moran_non_residualized)

# Plot 
ri_moran <- ggplot(df_moran_combined, aes(x = stat, fill = group)) +
  geom_histogram(binwidth = 0.00025, color = "darkslategrey", alpha = 0.6, position = "identity", show.legend = FALSE) +  
  geom_vline(aes(xintercept = obs_stat_moran_res, color = "Residualized", linetype = "Residualized"), size = .8) +
  geom_vline(aes(xintercept = obs_stat_moran_nr, color = "Non-residualized", linetype = "Non-residualized"), size = .8) +
  scale_color_manual(name = "Observed:", 
                     values = c("Residualized" = "black", "Non-residualized" = "deepskyblue4")) +
  scale_linetype_manual(name = "Observed:",
                        values = c("Residualized" = "solid", "Non-residualized" = "dashed")) +
  labs(title = "RI for Moran's I",
       x = NULL,
       y = "Frequency") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(),  # Removes major grid lines
        panel.grid.minor = element_blank(),
        legend.position = "bottom",
        legend.box = "horizontal",
        plot.title = element_text(size = 15), 
        axis.title.x = element_text(size = 14),  
        axis.title.y = element_text(size = 14),  
        axis.text = element_text(size = 13),     
        legend.text = element_text(size = 13),   
        legend.title = element_text(size = 14)) 

print(ri_moran)
ggsave("images/ri_moran.pdf",ri_moran, width = 10, height = 5)


# p-values RI, non-residualized vs. residualized
calc_two_tailed_p_value <- function(obs_stat_runs_res, randomized_stats) {
  extreme_count <- sum(abs(randomized_stats) >= abs(obs_stat_runs_res))   # Count how many random statistics are more extreme than the observed statistic
  p_value <- (extreme_count + 1) / (length(randomized_stats) + 1) # +1 for continuity correction
  return(p_value)
}

p_value_moran_nr_two_tailed <- calc_two_tailed_p_value(obs_stat_moran_nr, randomized_moran_non_residualized)
print(p_value_moran_nr_two_tailed)

p_value_moran_res_two_tailed <- calc_two_tailed_p_value(obs_stat_moran_res, randomized_moran_residualized)
print(p_value_moran_res_two_tailed)


################################################################################################################################
### Table A5: Understanding relationship between Delta diversity, Expellee share, Religious fractionalization in 1939 and 1950 #
################################################################################################################################

#Get rho of unconditional correlation mentioned on p. 11-12:
cor.test(dat$DeltaFrac, dat$srefugees1950) #0.25 (p<0.001)


f1<-"lm(DeltaFrac~srefugees1950, data=dat)"
f2<-"lm(DeltaFrac~srefugees1950+as.factor(kreis_kenn1950_alt), data=dat)"
f3<-"lm(DeltaFrac~srefugees1950+relfrac1939+sagri1939+sfemale1939+lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939, data=dat)"
f4<-"lm(DeltaFrac~srefugees1950+relfrac1939+sagri1939+sfemale1939+lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)"

f5<-"lm(relfrac1939~srefugees1950, data=dat)"
f6<-"lm(relfrac1939~srefugees1950+as.factor(kreis_kenn1950_alt), data=dat)"
f7<-"lm(relfrac1950~srefugees1950, data=dat)"
f8<-"lm(relfrac1950~srefugees1950 +as.factor(kreis_kenn1950_alt), data=dat)"


reg.formulas<-list(f1,f2,f3,f4,f5,f6, f7, f8)

models <- run_models(reg.formulas) 
models.cl <- run_coef(models)

TableA5panelA <- texreg(
  l = list(models[[1]], models[[2]], models[[3]],  models[[4]]), 
  override.se=list(models.cl[[1]][,2], models.cl[[2]][,2],  models.cl[[3]][,2], models.cl[[4]][,2]), 
  override.pval=list(models.cl[[1]][,4], models.cl[[2]][,4],  models.cl[[3]][,4], models.cl[[4]][,4]),     
  custom.coef.names=c("Share expellees (1950)", "Religious fractionalization (1939)","Share in agriculture (1939)", "Share female (1939)",  
                      "Ln distance to Eastern Border", "Distance to station (1950)", "Population density (1939)", "Destruction 1945", "Ln population (1939)"),
  custom.model.names = c("uncond", "FE", "Covars", "FE"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:deltafrac",
  omit.coef = "kreis_kenn1950_alt|Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE
)

TableA5panelB <- texreg(
  l = list(models[[5]], models[[6]], models[[7]], models[[8]]), 
  override.se=list(models.cl[[5]][,2], models.cl[[6]][,2],  models.cl[[7]][,2], models.cl[[8]][,2]), 
  override.pval=list(models.cl[[5]][,4], models.cl[[6]][,4],  models.cl[[7]][,4], models.cl[[8]][,4]),     
  custom.coef.names=c("Share expellees (1950)"),
  custom.model.names = c("uncond", "FE", "uncond", "FE"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:deltafrac",
  omit.coef = "kreis_kenn1950_alt|Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE
)

#########################################################################
###### Table 1: ATTENDANCE AT CATHOLIC MASS.  ###########################
#########################################################################


relig_dv<-list("lm(attendants_rel1970s~", "lm(attendants_rel1980s~", "lm(attendants_rel1990s~") 

expvars1 <- list("DeltaFrac+relfrac1939, data=dat)")

expvars2 <- list("DeltaFrac+ relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939, data=dat)")

expvars3 <- list("DeltaFrac+ relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)")

  
reg.formula1 <- formulas(relig_dv, expvars1)
reg.formula2 <- formulas(relig_dv, expvars2)
reg.formula3 <- formulas(relig_dv, expvars3)

mod1 <- run_models(c(reg.formula1, reg.formula2, reg.formula3)) 
mod.cl1 <- run_coef(mod1)


#DV means
means_dv <- c(rep(mean(dat$attendants_rel1970s,na.rm=T),3),rep(mean(dat$attendants_rel1980s,na.rm=T),3),rep(mean(dat$attendants_rel1990s,na.rm=T),3))
sds_dv <- c(rep(sd(dat$attendants_rel1970s,na.rm=T),3),rep(sd(dat$attendants_rel1980s,na.rm=T),3),rep(sd(dat$attendants_rel1990s,na.rm=T),3))

  
table1 <- texreg(
  l = list(mod1[[1]], mod1[[4]], mod1[[7]],  mod1[[2]], mod1[[5]], mod1[[8]], mod1[[3]], mod1[[6]], mod1[[9]]),
  override.se=list(mod.cl1[[1]][,2], mod.cl1[[4]][,2],  mod.cl1[[7]][,2], mod.cl1[[2]][,2],  mod.cl1[[5]][,2], mod.cl1[[8]][,2], mod.cl1[[3]][,2],  mod.cl1[[6]][,2], mod.cl1[[9]][,2]),
  override.pval=list(mod.cl1[[1]][,4], mod.cl1[[4]][,4],  mod.cl1[[7]][,4], mod.cl1[[2]][,4], mod.cl1[[5]][,4],mod.cl1[[8]][,4], mod.cl1[[3]][,4], mod.cl1[[6]][,4],mod.cl1[[9]][,4]),     
  custom.coef.names=c("Delta Diversity", "Religious fractionalization (1939)", "Share expellees","Share in agriculture (1939)", "Share female (1939)",  
                     "Ln distance to Eastern border", "Distance to station (1950)", "Population density (1939)", "Destruction 1945", "Ln population (1939)"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:selfemp",
  omit.coef = "kreis_kenn1950_alt|Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE, 
  custom.gof.rows = list("DV Mean" = means_dv, "DV SD"=sds_dv)
)

###################################################################################
###### Table A7 Panel A: CONLEY ERRORS for same specifications as in Table 1 ######
###################################################################################

#setting up formulas
outcomes = list("conleyreg(attendants_rel1970s~", "conleyreg(attendants_rel1980s~", "conleyreg(attendants_rel1990s~")
expvars1 = list("DeltaFrac+ relfrac1939, data=dat, 250, lat = 'lat', lon = 'lon')",
                "DeltaFrac+ relfrac1939+srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939, data=dat, 250, lat = 'lat', lon = 'lon')",
                "DeltaFrac+ relfrac1939+srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+county_id_f, data=dat, 250, lat = 'lat', lon = 'lon')")


#regressions
fs <- formulas(outcomes, expvars1)
ms <- run_models(fs)

#This table only changes standard errors, R2 and number of observations are from Table 1. 
texreg(
  l = ms,
  override.se=  list(ms[[1]][,2], ms[[2]][,2],  ms[[3]][,2], ms[[4]][,2],  ms[[5]][,2], ms[[6]][,2], ms[[7]][,2],  ms[[8]][,2], ms[[9]][,2]),
  override.pval=list(ms[[1]][,4], ms[[2]][,4],  ms[[3]][,4], ms[[4]][,4],  ms[[5]][,4], ms[[6]][,4], ms[[7]][,4],  ms[[8]][,4], ms[[9]][,4]),
  custom.coef.names = c('Delta diversity', "Religious fractionalization (1939)", "Share expellees",
                        "Share in agriculture (1939)", "Share female (1939)", "Ln distance to Eastern border", 
                        "Distance to station (1950)", "Population density (1939)", "Destruction", "Ln population (1939)"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  omit.coef = "county|Intercept",
)

###############################################################################
###### Table A7 Panel B: Alternative specification  w/out srefugees1950  ######
###############################################################################

relig_dv<-list("lm(attendants_rel1970s~", "lm(attendants_rel1980s~", "lm(attendants_rel1990s~") 

expvars1 <- list("DeltaFrac+ relfrac1939, data=dat)")

expvars2 <- list("DeltaFrac+ relfrac1939+ sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939, data=dat)")

expvars3 <- list("DeltaFrac+ relfrac1939+ sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)")


reg.formula1 <- formulas(relig_dv, expvars1)
reg.formula2 <- formulas(relig_dv, expvars2)
reg.formula3 <- formulas(relig_dv, expvars3)

mod1 <- run_models(c(reg.formula1, reg.formula2, reg.formula3)) 
mod.cl1 <- run_coef(mod1)


tableA7 <- texreg(
  l = list(mod1[[1]], mod1[[4]], mod1[[7]],  mod1[[2]], mod1[[5]], mod1[[8]], mod1[[3]], mod1[[6]], mod1[[9]]),
  override.se=list(mod.cl1[[1]][,2], mod.cl1[[4]][,2],  mod.cl1[[7]][,2], mod.cl1[[2]][,2],  mod.cl1[[5]][,2], mod.cl1[[8]][,2], mod.cl1[[3]][,2],  mod.cl1[[6]][,2], mod.cl1[[9]][,2]),
  override.pval=list(mod.cl1[[1]][,4], mod.cl1[[4]][,4],  mod.cl1[[7]][,4], mod.cl1[[2]][,4], mod.cl1[[5]][,4],mod.cl1[[8]][,4], mod.cl1[[3]][,4], mod.cl1[[6]][,4],mod.cl1[[9]][,4]),     
  custom.coef.names=c("Delta diversity", "Religious fractionalization (1939)", "Share in agriculture (1939)", "Share female (1939)",  
                     "Ln distance to Eastern border", "Distance to station (1950)", "Population density (1939)", "Destruction 1945", "Ln population (1939)"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:church_altspecs",
  omit.coef = "kreis_kenn1950_alt|Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE, 
)


########################################################################
######   Table A8:   MEMBERSHIP IN THE CHURCH: main specification  #####
########################################################################

church_dv<-list("lm(preligious_70s~", "lm(preligious_87s~") # 

expvars1 <- list("DeltaFrac+relfrac1939, data=dat)")

expvars2 <- list("DeltaFrac+ relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939, data=dat)")

expvars3 <- list("DeltaFrac+ relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)")

expvars4 <- list("DeltaFrac+ relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(diocese), data=dat)")


reg.formula1 <- formulas(church_dv, expvars1)
reg.formula2 <- formulas(church_dv, expvars2)
reg.formula3 <- formulas(church_dv, expvars3)
reg.formula4 <- formulas(church_dv, expvars4)

mod1<- run_models(c(reg.formula1, reg.formula2, reg.formula3, reg.formula4)) 
mod.cl1 <- run_coef(mod1)

#DV means
means_dv <- c(rep(mean(dat$preligious_70s,na.rm=T),4),rep(mean(dat$preligious_87s,na.rm=T),4))
sds_dv <- c(rep(sd(dat$preligious_70s,na.rm=T),4),rep(sd(dat$preligious_87s,na.rm=T),4))

church_affil_main <- texreg(
  l = list(mod1[[1]],  mod1[[3]], mod1[[5]], mod1[[7]], mod1[[2]],  mod1[[4]], mod1[[6]],  mod1[[8]]),
  override.se=list( mod.cl1[[1]][,2], mod.cl1[[3]][,2], mod.cl1[[5]][,2], mod.cl1[[7]][,2], mod.cl1[[2]][,2], mod.cl1[[4]][,2], mod.cl1[[6]][,2], mod.cl1[[8]][,2]),
  override.pval=list(mod.cl1[[1]][,4], mod.cl1[[3]][,4], mod.cl1[[5]][,4],mod.cl1[[7]][,4], mod.cl1[[2]][,4], mod.cl1[[4]][,4], mod.cl1[[6]][,4], mod.cl1[[8]][,4]),     
  custom.coef.names=c("Delta diversity", "Religious fractionalization (1939)","Share expellees (1950)", "Share in agriculture (1939)", "Share female (1939)",  
                      "Ln distance to Eastern Border", "Distance to station (1950)", "Population density (1939)", "Destruction", "Ln population (1939)"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:affiliationMain",
  omit.coef = "kreis_kenn1950_alt|diocese|Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE, 
  custom.gof.rows = list("DV Mean" = means_dv, "DV SD" = sds_dv)
)


#############################################################################################################
######   Table A9:   MEMBERSHIP IN THE CHURCH: alternative specification w/out srefugees1950  ####
#############################################################################################################
church_dv<-list("lm(preligious_70s~", "lm(preligious_87s~") # 

expvars1 <- list("DeltaFrac+relfrac1939, data=dat)")

expvars2 <- list("DeltaFrac+ relfrac1939+ sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939, data=dat)")

expvars3 <- list("DeltaFrac+ relfrac1939+ sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)")

expvars4 <- list("DeltaFrac+ relfrac1939+ sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(diocese), data=dat)")



reg.formula1 <- formulas(church_dv, expvars1)
reg.formula2 <- formulas(church_dv, expvars2)
reg.formula3 <- formulas(church_dv, expvars3)
reg.formula4 <- formulas(church_dv, expvars4)

mod1<- run_models(c(reg.formula1, reg.formula2, reg.formula3, reg.formula4)) 
mod.cl1 <- run_coef(mod1)

#DV means
means_dv <- c(rep(mean(dat$preligious_70s,na.rm=T),4),rep(mean(dat$preligious_87s,na.rm=T),4))
sds_dv <- c(rep(sd(dat$preligious_70s,na.rm=T),4),rep(sd(dat$preligious_87s,na.rm=T),4))

#Church membership
tableA9affil <- texreg(
  l = list(mod1[[1]],  mod1[[3]], mod1[[5]], mod1[[7]], mod1[[2]],  mod1[[4]], mod1[[6]],  mod1[[8]]),
  override.se=list( mod.cl1[[1]][,2], mod.cl1[[3]][,2], mod.cl1[[5]][,2], mod.cl1[[7]][,2], mod.cl1[[2]][,2], mod.cl1[[4]][,2], mod.cl1[[6]][,2], mod.cl1[[8]][,2]),
  override.pval=list(mod.cl1[[1]][,4], mod.cl1[[3]][,4], mod.cl1[[5]][,4],mod.cl1[[7]][,4], mod.cl1[[2]][,4], mod.cl1[[4]][,4], mod.cl1[[6]][,4], mod.cl1[[8]][,4]),     
  custom.coef.names=c("Delta diversity", "Religious fractionalization (1939)","Share in agriculture (1939)", "Share female (1939)",  
                      "Ln distance to Eastern border", "Distance to station (1950)", "Population density (1939)", "Destruction", "Ln population (1939)"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:affiliation",
  omit.coef = "kreis_kenn1950_alt|diocese|Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE, 
  custom.gof.rows = list("DV Mean" = means_dv, "DV SD" = sds_dv)
)


#########################################################################################################
###### Table A10 and Figure 3: CHURCH MEMBERSHIP IN PROTESTANT/CATHOLIC MUNICIPALITIES        ############
#########################################################################################################

protestant_municipalities<-subset(dat, sprotestants1939>0.9) #N=222
catholic_municipalities<-subset(dat, scatholics1939>0.9) 

church_dv<-list("lm(preligious_70s~", "lm(preligious_87s~")

expvars1 <- list("DeltaFrac+relfrac1939+srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939, data=catholic_municipalities)") 

expvars2 <- list("DeltaFrac+relfrac1939+srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=catholic_municipalities)")
expvars3 <- list("DeltaFrac+relfrac1939+srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939, data=protestant_municipalities)")

expvars4 <- list("DeltaFrac+relfrac1939+srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=protestant_municipalities)")


reg.formula1 <- formulas(church_dv, expvars1)
reg.formula2 <- formulas(church_dv, expvars2)
reg.formula3 <- formulas(church_dv, expvars3)
reg.formula4 <- formulas(church_dv, expvars4)

mod_affil<- run_models(c(reg.formula1, reg.formula2, reg.formula3, reg.formula4))
mod_affil.cl <- run_coef(mod_affil)


tableA10 <- texreg(
  l = list(mod_affil[[1]], mod_affil[[3]], mod_affil[[2]],  mod_affil[[4]], mod_affil[[5]], mod_affil[[7]], mod_affil[[6]], mod_affil[[8]]),
  override.se=list(mod_affil.cl[[1]][,2], mod_affil.cl[[3]][,2],  mod_affil.cl[[2]][,2], mod_affil.cl[[4]][,2], mod_affil.cl[[5]][,2], mod_affil.cl[[7]][,2], mod_affil.cl[[6]][,2], mod_affil.cl[[8]][,2]),
  override.pval=list(mod_affil.cl[[1]][,4], mod_affil.cl[[3]][,4],  mod_affil.cl[[2]][,4], mod_affil.cl[[4]][,4], mod_affil.cl[[5]][,4], mod_affil.cl[[7]][,4], mod_affil.cl[[6]][,4], mod_affil.cl[[8]][,4]),     
  custom.coef.names=c("Delta diversity", "Religious fractionalization (1939)","Share expellees", "Share in agriculture (1939)", "Share female (1939)",  
                      "Ln distance to Eastern border", "Distance to station (1950)", "Population density (1939)", "Destruction", "Ln population (1939)"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:affiliationProtCatholic", 
  omit.coef = "kreis_kenn1950_alt|diocese|Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE
)

#plotting this: 

plot.cath <- sapply(list(mod_affil[[3]], mod_affil[[4]]), function(x){
  coeftest(x, vcov.=hc1)[2,1:2] 
})

plot.prot <- sapply(list(mod_affil[[7]], mod_affil[[8]]), function(x){
  coeftest(x, vcov.=hc1)[2,1:2] 
})


plot.out <- as.data.frame(t(plot.cath))
plot.out$DV<-c("1970", "1987")
plot.out$IV<-rep("Delta Diversity", 2)
plot.out$lb <- plot.out$Estimate - 1.96*plot.out$`Std. Error`
plot.out$ub <- plot.out$Estimate + 1.96*plot.out$`Std. Error`


##Figure 3 left panel
figureAffiliationCath <- ggplot(plot.out, aes(x = DV, y = Estimate, shape = IV)) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_pointrange(aes(ymin = lb, ymax = ub), alpha = 0.75, position = position_dodge(width = 0.1), size = 1) +
  scale_y_continuous(name = "Average marginal effect: Church affiliation in Catholic municipalities") +
  labs(shape = "IV") +  # Change legend label to "IV"
  theme_bw() +
  theme(
    text = element_text(size = 12, family = "sans", color = "black"),
    legend.position = "bottom",
    legend.text = element_text(size = 12, family = "sans"),
    axis.text.x = element_text(angle = 47, hjust = 1, size = 12, family = "sans", color = "black"),
    axis.title.x = element_blank(),  # Remove x-axis label
    legend.title = element_blank(),
    plot.margin = unit(c(.5, .5, .5, 1), "cm")
  )

#ggsave("images/Figure3left.jpg", figureAffiliationCath, width = 8, height = 6)
ggsave("images/Figure3left.tiff", figureAffiliationCath, width = 8, height = 6, dpi = 300, device = "tiff")

plot.out2 <- as.data.frame(t(plot.prot))
plot.out2$DV<-c("1970", "1987")
plot.out2$IV<-rep("Delta Diversity", 2)
plot.out2$lb <- plot.out2$Estimate - 1.96*plot.out2$`Std. Error`
plot.out2$ub <- plot.out2$Estimate + 1.96*plot.out2$`Std. Error`

figureAffiliationProt <- ggplot(plot.out2, aes(x = DV, y = Estimate, shape = IV)) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_pointrange(aes(ymin = lb, ymax = ub), alpha = 0.75, position = position_dodge(width = 0.1), size = 1) +
  scale_y_continuous(name = "Average marginal effect: Church affiliation in Protestant municipalities") +
  labs(shape = "IV") +  # Change legend label to "IV"
  theme_bw() +
  theme(
    text = element_text(size = 12, family = "sans", color = "black"),
    legend.position = "bottom",
    legend.text = element_text(size = 12, family = "sans"),
    axis.text.x = element_text(angle = 47, hjust = 1, size = 12, family = "sans", color = "black"),
    axis.title.x = element_blank(),  # Remove x-axis label
    legend.title = element_blank(),
    plot.margin = unit(c(.5, .5, .5, 1), "cm")
  )

#ggsave("images/Figure3right.jpg", figureAffiliationProt, width = 8, height = 6)

ggsave("images/Figure3right.tiff", figureAffiliationProt, width = 8, height = 6, dpi = 300, device = "tiff")

###########################################################################################################################
###### Table A11: CHURCH MEMBERSHIP IN PROTESTANT/CATHOLIC MUNICIPALITIES.  Alternative specification w/out srefugees1950
###########################################################################################################################


expvars1 <- list("DeltaFrac+relfrac1939+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939, data=catholic_municipalities)") 
expvars2 <- list("DeltaFrac+relfrac1939+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=catholic_municipalities)")
expvars3 <- list("DeltaFrac+relfrac1939+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939, data=protestant_municipalities)")
expvars4 <- list("DeltaFrac+relfrac1939+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=protestant_municipalities)")


reg.formula1 <- formulas(church_dv, expvars1)
reg.formula2 <- formulas(church_dv, expvars2)
reg.formula3 <- formulas(church_dv, expvars3)
reg.formula4 <- formulas(church_dv, expvars4)

mod_affil<- run_models(c(reg.formula1, reg.formula2, reg.formula3, reg.formula4)) 
mod_affil.cl <- run_coef(mod_affil)

#Church membership
tableA11 <- texreg(
  l = list(mod_affil[[1]], mod_affil[[3]], mod_affil[[2]],  mod_affil[[4]], mod_affil[[5]], mod_affil[[7]], mod_affil[[6]], mod_affil[[8]]),
  override.se=list(mod_affil.cl[[1]][,2], mod_affil.cl[[3]][,2],  mod_affil.cl[[2]][,2], mod_affil.cl[[4]][,2], mod_affil.cl[[5]][,2], mod_affil.cl[[7]][,2], mod_affil.cl[[6]][,2], mod_affil.cl[[8]][,2]),
  override.pval=list(mod_affil.cl[[1]][,4], mod_affil.cl[[3]][,4],  mod_affil.cl[[2]][,4], mod_affil.cl[[4]][,4], mod_affil.cl[[5]][,4], mod_affil.cl[[7]][,4], mod_affil.cl[[6]][,4], mod_affil.cl[[8]][,4]),     
  custom.coef.names=c("Delta diversity", "Religious fractionalization (1939)","Share in agriculture (1939)", "Share female (1939)",  
                      "Ln distance to Eastern border", "Distance to station (1950)", "Population density (1939)", "Destruction", "Ln population (1939)"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:affiliationProtCatholic_v2",
  omit.coef = "kreis_kenn1950_alt|diocese|Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE
)



####################################################################################
############## Table 2: DIVORCE RATES, Main specification ##########################
####################################################################################

divorce_dv<-list("lm(pdivorce_tomar_70~", "lm(pdivorce_tomar_87~")


expvars1 <- list("DeltaFrac+relfrac1939, data=dat)")

expvars2 <- list("DeltaFrac+ relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939, data=dat)")

expvars3 <- list("DeltaFrac+ relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)")

expvars4 <- list("DeltaFrac+ relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(diocese), data=dat)")


reg.formula1 <- formulas(divorce_dv, expvars1)
reg.formula2 <- formulas(divorce_dv, expvars2)
reg.formula3 <- formulas(divorce_dv, expvars3)
reg.formula4 <- formulas(divorce_dv, expvars4)

mod_divorce <- run_models(c(reg.formula1, reg.formula2, reg.formula3, reg.formula4)) 
mod.divorcecl1 <- run_coef(mod_divorce)

##DV means
means_dv <- c(rep(mean(dat$pdivorce_tomar_70,na.rm=T),4),rep(mean(dat$pdivorce_tomar_87,na.rm=T), 4))
sds_dv <- c(rep(sd(dat$pdivorce_tomar_70,na.rm=T),4),rep(sd(dat$pdivorce_tomar_87,na.rm=T),4))


table2 <- texreg(
  l = list(mod_divorce[[1]],  mod_divorce[[3]], mod_divorce[[5]], mod_divorce[[7]], mod_divorce[[2]],  mod_divorce[[4]], mod_divorce[[6]],  mod_divorce[[8]]),
  override.se=list(mod.divorcecl1[[1]][,2], mod.divorcecl1[[3]][,2], mod.divorcecl1[[5]][,2], mod.divorcecl1[[7]][,2], mod.divorcecl1[[2]][,2], mod.divorcecl1[[4]][,2], mod.divorcecl1[[6]][,2], mod.divorcecl1[[8]][,2]),
  override.pval=list(mod.divorcecl1[[1]][,4], mod.divorcecl1[[3]][,4], mod.divorcecl1[[5]][,4],mod.divorcecl1[[7]][,4], mod.divorcecl1[[2]][,4], mod.divorcecl1[[4]][,4], mod.divorcecl1[[6]][,4], mod.divorcecl1[[8]][,4]),     
  custom.coef.names=c("Delta diversity", "Religious fractionalization (1939)", "Share expellees", "Share in agriculture (1939)", "Share female (1939)",  
                      "Ln distance to Eastern border", "Distance to station (1950)", "Population density (1939)","Destruction", "Ln population (1939)"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:divorce",
  omit.coef = "kreis_kenn1950_alt|diocese|Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE, 
  custom.gof.rows = list("DV Mean" = means_dv, "DV SD" = sds_dv)
)



################################################################################################
###### Table A12, Panel A: Religious diversity and divorce rates with Conley errors ############
################################################################################################

divorce_dv<-list("conleyreg(pdivorce_tomar_70~", "conleyreg(pdivorce_tomar_87~")
expvars1 = list("DeltaFrac+ relfrac1939, data=dat, 250, lat = 'lat', lon = 'lon')",
                "DeltaFrac+ relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939, data=dat, 250, lat = 'lat', lon = 'lon')",
                "DeltaFrac+ relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+county_id_f, data=dat, 250, lat = 'lat', lon = 'lon')",
                "DeltaFrac+ relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+diocese_id_f, data=dat, 250, lat = 'lat', lon = 'lon')")

#regressions
fs <- formulas(divorce_dv, expvars1)
ms <- run_models(fs)


#Number of observations and adjusted R2 in this table is from Table 2
texreg(
  l = ms,
  override.se=  list(ms[[1]][,2], ms[[2]][,2],  ms[[3]][,2], ms[[4]][,2],  ms[[5]][,2], ms[[6]][,2], ms[[7]][,2], ms[[8]][,2]),
  override.pval=list(ms[[1]][,4], ms[[2]][,4],  ms[[3]][,4], ms[[4]][,4],  ms[[5]][,4], ms[[6]][,4], ms[[7]][,4], ms[[8]][,4]),
 custom.coef.names = c("Delta diversity", "Religious fractionalization (1939)","Share expellees",
  "Share in agriculture (1939)", "Share female (1939)", "Ln distance to Eastern border", 
  "Distance to station (1950)", "Population density (1939)","Destruction", "Ln population (1939)"), 
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  omit.coef = "diocese|county|Intercept")

##########################################################################################################
############## Table A12, Panel B: DIVORCE RATES: Alternative specification without expellee share #########
##########################################################################################################

divorce_dv<-list("lm(pdivorce_tomar_70~", "lm(pdivorce_tomar_87~")


expvars1 <- list("DeltaFrac+relfrac1939, data=dat)")

expvars2 <- list("DeltaFrac+ relfrac1939+ sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939, data=dat)")

expvars3 <- list("DeltaFrac+ relfrac1939+ sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)")

expvars4 <- list("DeltaFrac+ relfrac1939+ sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(diocese), data=dat)")


reg.formula1 <- formulas(divorce_dv, expvars1)
reg.formula2 <- formulas(divorce_dv, expvars2)
reg.formula3 <- formulas(divorce_dv, expvars3)
reg.formula4 <- formulas(divorce_dv, expvars4)

mod_divorce <- run_models(c(reg.formula1, reg.formula2, reg.formula3, reg.formula4)) 
mod.divorcecl1 <- run_coef(mod_divorce)


tableA12b <- texreg(
  l = list(mod_divorce[[1]],  mod_divorce[[3]], mod_divorce[[5]], mod_divorce[[7]], mod_divorce[[2]],  mod_divorce[[4]], mod_divorce[[6]],  mod_divorce[[8]]),
  override.se=list(mod.divorcecl1[[1]][,2], mod.divorcecl1[[3]][,2], mod.divorcecl1[[5]][,2], mod.divorcecl1[[7]][,2], mod.divorcecl1[[2]][,2], mod.divorcecl1[[4]][,2], mod.divorcecl1[[6]][,2], mod.divorcecl1[[8]][,2]),
  override.pval=list(mod.divorcecl1[[1]][,4], mod.divorcecl1[[3]][,4], mod.divorcecl1[[5]][,4],mod.divorcecl1[[7]][,4], mod.divorcecl1[[2]][,4], mod.divorcecl1[[4]][,4], mod.divorcecl1[[6]][,4], mod.divorcecl1[[8]][,4]),     
  custom.coef.names=c("Delta diversity", "Religious fractionalization (1939)",  "Share in agriculture (1939)", "Share female (1939)",  
                      "Ln distance to Eastern border", "Distance to station (1950)", "Population density (1939)","Destruction", "Ln population (1939)"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:divorceAlt",
  omit.coef = "kreis_kenn1950_alt|diocese|Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE, 
)



#########################################################################################################
################ Table A13, panel A: Divorce rates in subsets of Catholic & Protestant Municipalities ###
#########################################################################################################



divorce_dv<-list("lm(pdivorce_tomar_70~", "lm(pdivorce_tomar_87~")

expvars1 <- list("DeltaFrac+relfrac1939+sagri1939+srefugees1950+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939, data=catholic_municipalities)")

expvars2 <- list("DeltaFrac+relfrac1939+sagri1939+srefugees1950+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=catholic_municipalities)")

expvars3 <- list("DeltaFrac+relfrac1939+sagri1939+srefugees1950+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939, data=protestant_municipalities)")

expvars4 <- list("DeltaFrac+relfrac1939+sagri1939+srefugees1950+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=protestant_municipalities)")

reg.formula1 <- formulas(divorce_dv, expvars1)
reg.formula2 <- formulas(divorce_dv, expvars2)
reg.formula3 <- formulas(divorce_dv, expvars3)
reg.formula4 <- formulas(divorce_dv, expvars4)

mod_divorce<- run_models(c(reg.formula1, reg.formula2, reg.formula3, reg.formula4)) 
mod_divorce.cl <- run_coef(mod_divorce)

##DV means
means_dv <- c(rep(mean(catholic_municipalities$pdivorce_tomar_70,na.rm=T),2),rep(mean(catholic_municipalities$pdivorce_tomar_87,na.rm=T), 2), rep(mean(protestant_municipalities$pdivorce_tomar_70,na.rm=T),2),rep(mean(protestant_municipalities$pdivorce_tomar_87,na.rm=T), 2))
sds_dv <- c(rep(sd(catholic_municipalities$pdivorce_tomar_70,na.rm=T),2),rep(sd(catholic_municipalities$pdivorce_tomar_87,na.rm=T),2),rep(sd(protestant_municipalities$pdivorce_tomar_70,na.rm=T),2), rep(sd(protestant_municipalities$pdivorce_tomar_87,na.rm=T), 2))


#Divorce rate
tableA13a <- texreg(
  l = list(mod_divorce[[1]], mod_divorce[[3]], mod_divorce[[2]],  mod_divorce[[4]], mod_divorce[[5]], mod_divorce[[7]], mod_divorce[[6]], mod_divorce[[8]]),
  override.se=list(mod_divorce.cl[[1]][,2], mod_divorce.cl[[3]][,2],  mod_divorce.cl[[2]][,2], mod_divorce.cl[[4]][,2], mod_divorce.cl[[5]][,2], mod_divorce.cl[[7]][,2], mod_divorce.cl[[6]][,2], mod_divorce.cl[[8]][,2]),
  override.pval=list(mod_divorce.cl[[1]][,4], mod_divorce.cl[[3]][,4],  mod_divorce.cl[[2]][,4], mod_divorce.cl[[4]][,4], mod_divorce.cl[[5]][,4], mod_divorce.cl[[7]][,4], mod_divorce.cl[[6]][,4], mod_divorce.cl[[8]][,4]),     
  custom.coef.names=c("Delta diversity", "Religious fractionalization (1939)", "Share expellees","Share in agriculture (1939)", "Share female (1939)",  
                      "Ln distance to Eastern border", "Distance to station (1950)",  "Population density (1939)", "Destruction", "Ln population (1939)"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:divorces-separate",
  omit.coef = "kreis_kenn1950_alt|diocese|Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE
)

#########################
#### Creating Figure 4: 
#########################

plot.cath <- sapply(list(mod_divorce[[3]], mod_divorce[[4]]), function(x){
  coeftest(x, vcov.=hc1)[2,1:2] 
})

plot.prot <- sapply(list(mod_divorce[[7]], mod_divorce[[8]]), function(x){
  coeftest(x, vcov.=hc1)[2,1:2] 
})


plot.out <- as.data.frame(t(plot.cath))
plot.out$DV<-c("1970", "1987")
plot.out$IV<-rep("Delta diversity", 2)
plot.out$lb <- plot.out$Estimate - 1.96*plot.out$`Std. Error`
plot.out$ub <- plot.out$Estimate + 1.96*plot.out$`Std. Error`


plot.out2 <- as.data.frame(t(plot.prot))
plot.out2$DV<-c("1970", "1987")
plot.out2$IV<-rep("Delta diversity", 2)
plot.out2$lb <- plot.out2$Estimate - 1.96*plot.out2$`Std. Error`
plot.out2$ub <- plot.out2$Estimate + 1.96*plot.out2$`Std. Error`


figureDivorceCath <- ggplot(plot.out, aes(x = DV, y = Estimate, shape = IV)) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_pointrange(aes(ymin = lb, ymax = ub), alpha = 0.75, position = position_dodge(width = 0.1), size = 1) +
  scale_y_continuous(name = "Average marginal effect: Divorce rates in Catholic municipalities") +
  labs(shape = "IV") +  # Change legend label to "IV"
  theme_bw() +
  theme(
    text = element_text(size = 12, family = "sans", color = "black"),
    legend.position = "bottom",
    legend.text = element_text(size = 12, family = "sans"),
    axis.text.x = element_text(angle = 47, hjust = 1, size = 12, family = "sans", color = "black"),
    axis.title.x = element_blank(),  # Remove x-axis label
    legend.title = element_blank(),
    plot.margin = unit(c(.5, .5, .5, 1), "cm")
  )

ggsave("images/Figure4left.tiff", figureDivorceCath, width = 8, height = 6, dpi = 300, device = "tiff")

figureDivorceProt <- ggplot(plot.out2, aes(x = DV, y = Estimate, shape = IV)) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_pointrange(aes(ymin = lb, ymax = ub), alpha = 0.75, position = position_dodge(width = 0.1), size = 1) +
  scale_y_continuous(name = "Average marginal effect: Divorce rates in Protestant municipalities") +
  labs(shape = "IV") +  # Change legend label to "IV"
  theme_bw() +
  theme(
    text = element_text(size = 12, family = "sans", color = "black"),
    legend.position = "bottom",
    legend.text = element_text(size = 12, family = "sans"),
    axis.text.x = element_text(angle = 47, hjust = 1, size = 12, family = "sans", color = "black"),
    axis.title.x = element_blank(),  # Remove x-axis label
    legend.title = element_blank(),
    plot.margin = unit(c(.5, .5, .5, 1), "cm")
  )

ggsave("images/Figure4right.tiff", figureDivorceProt, width = 8, height = 6, dpi = 300, device = "tiff")

#################################################################################################################################
################ Table A13, Panel B: Divorce rates in subsets of Catholic & Protestant Municipalities without expellee share ####
#################################################################################################################################

divorce_dv<-list("lm(pdivorce_tomar_70~", "lm(pdivorce_tomar_87~")

expvars1 <- list("DeltaFrac+relfrac1939+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939, data=catholic_municipalities)")

expvars2 <- list("DeltaFrac+relfrac1939+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=catholic_municipalities)")

expvars3 <- list("DeltaFrac+relfrac1939+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939, data=protestant_municipalities)")

expvars4 <- list("DeltaFrac+relfrac1939+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=protestant_municipalities)")

reg.formula1 <- formulas(divorce_dv, expvars1)
reg.formula2 <- formulas(divorce_dv, expvars2)
reg.formula3 <- formulas(divorce_dv, expvars3)
reg.formula4 <- formulas(divorce_dv, expvars4)

mod_divorce<- run_models(c(reg.formula1, reg.formula2, reg.formula3, reg.formula4)) 
mod_divorce.cl <- run_coef(mod_divorce)

##DV means
means_dv <- c(rep(mean(catholic_municipalities$pdivorce_tomar_70,na.rm=T),2),rep(mean(catholic_municipalities$pdivorce_tomar_87,na.rm=T), 2), rep(mean(protestant_municipalities$pdivorce_tomar_70,na.rm=T),2),rep(mean(protestant_municipalities$pdivorce_tomar_87,na.rm=T), 2))
sds_dv <- c(rep(sd(catholic_municipalities$pdivorce_tomar_70,na.rm=T),2),rep(sd(catholic_municipalities$pdivorce_tomar_87,na.rm=T),2),rep(sd(protestant_municipalities$pdivorce_tomar_70,na.rm=T),2), rep(sd(protestant_municipalities$pdivorce_tomar_87,na.rm=T), 2))


#Divorce rate
tableA13b <- texreg(
  l = list(mod_divorce[[1]], mod_divorce[[3]], mod_divorce[[2]],  mod_divorce[[4]], mod_divorce[[5]], mod_divorce[[7]], mod_divorce[[6]], mod_divorce[[8]]),
  override.se=list(mod_divorce.cl[[1]][,2], mod_divorce.cl[[3]][,2],  mod_divorce.cl[[2]][,2], mod_divorce.cl[[4]][,2], mod_divorce.cl[[5]][,2], mod_divorce.cl[[7]][,2], mod_divorce.cl[[6]][,2], mod_divorce.cl[[8]][,2]),
  override.pval=list(mod_divorce.cl[[1]][,4], mod_divorce.cl[[3]][,4],  mod_divorce.cl[[2]][,4], mod_divorce.cl[[4]][,4], mod_divorce.cl[[5]][,4], mod_divorce.cl[[7]][,4], mod_divorce.cl[[6]][,4], mod_divorce.cl[[8]][,4]),     
  custom.coef.names=c("Delta Diversity", "Religious Fractionalization (1939)","Share in agriculture (1939)", "Share female (1939)",  
                      "Ln Distance to Eastern Border", "Distance to station (1950)",  "Population density (1939)", "Destruction", "Ln population (1939)"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:divorceAlt",
  omit.coef = "kreis_kenn1950_alt|diocese|Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE, 
  custom.gof.rows = list("DV Mean" = means_dv, "DV SD" = sds_dv)
)


###################################################################################################
######## Table A15: Support for CSU  1953-1990  with refugee control and county fixed effects #####
###################################################################################################


dvCSU<-list("lm(csu_voters1953~", "lm(csu_voters1957~", "lm(csu_voters1961~", "lm(csu_voters1965~", "lm(csu_voters1969~", "lm(csu_voters1972~", 
            "lm(csu_voters1976~", "lm(csu_voters1980~", "lm(csu_voters1983~","lm(csu_voters1987~", "lm(csu_voters1990~")

expvars <- list("DeltaFrac+ relfrac1939+ srefugees1950+ sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)")

#create and run formulas + get coefficients
reg.formula2 <- formulas(dvCSU, expvars)
mod2 <- run_models(reg.formula2)
mod.cl2 <- run_coef(mod2)

##DV means
means_dv<-round(c(mean(dat$csu_voters1953, na.rm=TRUE), mean(dat$csu_voters1957, na.rm=TRUE), mean(dat$csu_voters1961, na.rm=TRUE), mean(dat$csu_voters1965, na.rm=TRUE), 
                  mean(dat$csu_voters1969, na.rm=TRUE), mean(dat$csu_voters1972, na.rm=TRUE), mean(dat$csu_voters1976, na.rm=TRUE), mean(dat$csu_voters1980, na.rm=TRUE), mean(dat$csu_voters1983, na.rm=TRUE), mean(dat$csu_voters1987, na.rm=TRUE), mean(dat$csu_voters1990, na.rm=TRUE)), 3)
sds_dv  <-round(c(sd(dat$csu_voters1953, na.rm=TRUE), sd(dat$csu_voters1957, na.rm=TRUE), sd(dat$csu_voters1961, na.rm=TRUE), sd(dat$csu_voters1965, na.rm=TRUE), 
                  sd(dat$csu_voters1969, na.rm=TRUE), sd(dat$csu_voters1972, na.rm=TRUE), sd(dat$csu_voters1976, na.rm=TRUE), sd(dat$csu_voters1980, na.rm=TRUE), sd(dat$csu_voters1983, na.rm=TRUE), sd(dat$csu_voters1987, na.rm=TRUE), sd(dat$csu_voters1990, na.rm=TRUE)), 3)

csu_table_main <- texreg(
  l = list(mod2[[1]], mod2[[2]], mod2[[3]], mod2[[4]], mod2[[5]], mod2[[6]], mod2[[7]], mod2[[8]], mod2[[9]], mod2[[10]], mod2[[11]]),
  override.se=list(mod.cl2[[1]][,2], mod.cl2[[2]][,2], mod.cl2[[3]][,2], mod.cl2[[4]][,2], mod.cl2[[5]][,2], mod.cl2[[6]][,2], 
                   mod.cl2[[7]][,2], mod.cl2[[8]][,2], mod.cl2[[9]][,2], mod.cl2[[10]][,2], mod.cl2[[11]][,2]),
  override.pval=list(mod.cl2[[1]][,4], mod.cl2[[2]][,4], mod.cl2[[3]][,4], mod.cl2[[4]][,4], mod.cl2[[5]][,4], mod.cl2[[6]][,4], 
                     mod.cl2[[7]][,4], mod.cl2[[8]][,4], mod.cl2[[9]][,4], mod.cl2[[10]][,4], mod.cl2[[11]][,4]),     
  custom.model.names = c("1953","1957","1961","1965","1969", "1972", "1976", "1980", "1983", "1987", "1990"),
  custom.coef.names=c("Delta diversity", "Religious fractionalization (1939)",  "Share expellees (1950)","Share in agriculture (1939)", "Share female (1939)",  
                      "Ln distance to Eastern Border", "Distance to station (1950)", "Population density (1939)","Destruction" ,"Ln population (1939)"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:CSUvote",
  omit.coef = "kreis_kenn1950_alt|Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE,  
  custom.gof.rows = list("DV Mean" = means_dv, "DV SD" = sds_dv)
)

#################################################
###############  Figure 5, left panel (CSU) #####
#################################################
plot.outCSURef <- sapply(list(mod2[[1]], mod2[[2]], mod2[[3]], mod2[[4]], mod2[[5]], mod2[[6]], mod2[[7]], mod2[[8]], mod2[[9]], mod2[[10]], mod2[[11]]), function(x){
  coeftest(x, vcov.=hc1)[2,1:2] 
})

plot.outCSURef <- as.data.frame(t(plot.outCSURef))
plot.outCSURef$DV<-c("1953","1957","1961","1965","1969", "1972", "1976", "1980", "1983", "1987", "1990")
plot.outCSURef$IV<-rep("Delta diversity", 11)
plot.outCSURef$lb <- plot.outCSURef$Estimate - 1.96*plot.outCSURef$`Std. Error`
plot.outCSURef$ub <- plot.outCSURef$Estimate + 1.96*plot.outCSURef$`Std. Error`


figureCSURef <- ggplot(plot.outCSURef, aes(x = DV, y = Estimate, color = IV)) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_pointrange(aes(ymin = lb, ymax = ub), 
                  alpha = 0.75, 
                  position = position_dodge(width = 0.1), 
                  size = 0.5) +  # Set size to 0.5 for smaller circles
  scale_y_continuous(name = "Average marginal effect: CSU vote") +
  labs(color = "IV") +
  scale_color_manual(values = "black") +  # Set symbol color to black
  theme_bw() +
  theme(
    text = element_text(size = 12, family = "sans", color = "black"),
    legend.position = "bottom",
    legend.text = element_text(size = 12, family = "sans"),
    axis.text.x = element_text(angle = 47, hjust = 1, size = 12, family = "sans", color = "black"),
    axis.title.x = element_blank(),  # Remove x-axis label
    legend.title = element_blank(),
    plot.margin = unit(c(.5, .5, .5, 1), "cm")
  )

ggsave("images/Figure5left.tiff",figureCSURef, width = 8, height = 6, dpi = 300, device = "tiff")

##############################################################################################
######## Table A16, panel A: Support for CSU  1953-1990: Alternative spec.  without FEs ######
##############################################################################################

dvCSU<-list("lm(csu_voters1953~", "lm(csu_voters1957~", "lm(csu_voters1961~", "lm(csu_voters1965~", "lm(csu_voters1969~", "lm(csu_voters1972~", 
            "lm(csu_voters1976~", "lm(csu_voters1980~", "lm(csu_voters1983~","lm(csu_voters1987~", "lm(csu_voters1990~")

expvars <- list("DeltaFrac+ relfrac1939+ srefugees1950+ sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939, data=dat)")

reg.formula2 <- formulas(dvCSU, expvars)
mod2 <- run_models(reg.formula2)
mod.cl2 <- run_coef(mod2)


csu_table_alt1 <- texreg(
  l = list(mod2[[1]], mod2[[2]], mod2[[3]], mod2[[4]], mod2[[5]], mod2[[6]], mod2[[7]], mod2[[8]], mod2[[9]], mod2[[10]], mod2[[11]]),
  override.se=list(mod.cl2[[1]][,2], mod.cl2[[2]][,2], mod.cl2[[3]][,2], mod.cl2[[4]][,2], mod.cl2[[5]][,2], mod.cl2[[6]][,2], 
                   mod.cl2[[7]][,2], mod.cl2[[8]][,2], mod.cl2[[9]][,2], mod.cl2[[10]][,2], mod.cl2[[11]][,2]),
  override.pval=list(mod.cl2[[1]][,4], mod.cl2[[2]][,4], mod.cl2[[3]][,4], mod.cl2[[4]][,4], mod.cl2[[5]][,4], mod.cl2[[6]][,4], 
                     mod.cl2[[7]][,4], mod.cl2[[8]][,4], mod.cl2[[9]][,4], mod.cl2[[10]][,4], mod.cl2[[11]][,4]),     
  custom.model.names = c("1953","1957","1961","1965","1969", "1972", "1976", "1980", "1983", "1987", "1990"),
  custom.coef.names=c("Delta diversity", "Religious fractionalization (1939)", "Share expellees", "Share in agriculture (1939)", "Share female (1939)",  
                      "Ln distance to Eastern Border", "Distance to station (1950)", "Population density (1939)","Destruction" ,"Ln population (1939)"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:CSUvoteAlt",
  omit.coef = "kreis_kenn1950_alt|Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE,   
)

#############################################################################################
######## Table A16, panel B: Support for CSU  1953-1990  with FEs for Catholic dioceses #####
#############################################################################################

dvCSU<-list("lm(csu_voters1953~", "lm(csu_voters1957~", "lm(csu_voters1961~", "lm(csu_voters1965~", "lm(csu_voters1969~", "lm(csu_voters1972~", 
            "lm(csu_voters1976~", "lm(csu_voters1980~", "lm(csu_voters1983~","lm(csu_voters1987~", "lm(csu_voters1990~")

expvars <- list("DeltaFrac+ relfrac1939+srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(diocese), data=dat)")

#create and run formulas + get coefficients
reg.formula2 <- formulas(dvCSU, expvars)
mod2 <- run_models(reg.formula2)
mod.cl2 <- run_coef(mod2)


csu_tableAlt2 <- texreg(
  l = list(mod2[[1]], mod2[[2]], mod2[[3]], mod2[[4]], mod2[[5]], mod2[[6]], mod2[[7]], mod2[[8]], mod2[[9]], mod2[[10]], mod2[[11]]),
  override.se=list(mod.cl2[[1]][,2], mod.cl2[[2]][,2], mod.cl2[[3]][,2], mod.cl2[[4]][,2], mod.cl2[[5]][,2], mod.cl2[[6]][,2], 
                   mod.cl2[[7]][,2], mod.cl2[[8]][,2], mod.cl2[[9]][,2], mod.cl2[[10]][,2], mod.cl2[[11]][,2]),
  override.pval=list(mod.cl2[[1]][,4], mod.cl2[[2]][,4], mod.cl2[[3]][,4], mod.cl2[[4]][,4], mod.cl2[[5]][,4], mod.cl2[[6]][,4], 
                     mod.cl2[[7]][,4], mod.cl2[[8]][,4], mod.cl2[[9]][,4], mod.cl2[[10]][,4], mod.cl2[[11]][,4]),     
  custom.model.names = c("1953","1957","1961","1965","1969", "1972", "1976", "1980", "1983", "1987", "1990"),
  custom.coef.names=c("Delta diversity", "Religious fractionalization (1939)", "Share expellees", "Share in agriculture (1939)", "Share female (1939)",  
                      "Ln distance to Eastern border", "Distance to station (1950)", "Population density (1939)","Destruction" ,"Ln population (1939)"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:CSUvoteAlt",
  omit.coef = "diocese|Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE,   
)


###################################################################################################################
######## Table A15, panel C: Support for CSU  1953-1990: Alternative specification 3, without srefugees1950 #######
###################################################################################################################


dvCSU<-list("lm(csu_voters1953~", "lm(csu_voters1957~", "lm(csu_voters1961~", "lm(csu_voters1965~", "lm(csu_voters1969~", "lm(csu_voters1972~", 
            "lm(csu_voters1976~", "lm(csu_voters1980~", "lm(csu_voters1983~","lm(csu_voters1987~", "lm(csu_voters1990~")

expvars <- list("DeltaFrac+ relfrac1939+ sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)")


reg.formula2 <- formulas(dvCSU, expvars)
mod2 <- run_models(reg.formula2)
mod.cl2 <- run_coef(mod2)


csu_table_Alt3 <- texreg(
  l = list(mod2[[1]], mod2[[2]], mod2[[3]], mod2[[4]], mod2[[5]], mod2[[6]], mod2[[7]], mod2[[8]], mod2[[9]], mod2[[10]], mod2[[11]]),
  override.se=list(mod.cl2[[1]][,2], mod.cl2[[2]][,2], mod.cl2[[3]][,2], mod.cl2[[4]][,2], mod.cl2[[5]][,2], mod.cl2[[6]][,2], 
                   mod.cl2[[7]][,2], mod.cl2[[8]][,2], mod.cl2[[9]][,2], mod.cl2[[10]][,2], mod.cl2[[11]][,2]),
  override.pval=list(mod.cl2[[1]][,4], mod.cl2[[2]][,4], mod.cl2[[3]][,4], mod.cl2[[4]][,4], mod.cl2[[5]][,4], mod.cl2[[6]][,4], 
                     mod.cl2[[7]][,4], mod.cl2[[8]][,4], mod.cl2[[9]][,4], mod.cl2[[10]][,4], mod.cl2[[11]][,4]),     
  custom.model.names = c("1953","1957","1961","1965","1969", "1972", "1976", "1980", "1983", "1987", "1990"),
  custom.coef.names=c("Delta diversity", "Religious fractionalization (1939)", "Share in agriculture (1939)", "Share female (1939)",
                     "Ln Distance to Eastern Border", "Distance to station (1950)", "Population density (1939)","Destruction" ,"Ln population (1939)"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:CSUvoteAlt",
  omit.coef = "kreis_kenn1950_alt|Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE)




###########################################################################
######## Table A17: Support for Greens  1980-1990: Main specification #####
###########################################################################

dvGreens<-list("lm(greens_voters1980~", "lm(greens_voters1983~", "lm(greens_voters1987~", "lm(greens_voters1990~")

expvars1 <- list("DeltaFrac+relfrac1939, data=dat)")

expvars2 <- list("DeltaFrac+ relfrac1939+srefugees1950+ sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)")

reg.formula1 <- formulas(dvGreens, expvars1)
reg.formula2 <- formulas(dvGreens, expvars2)

mod1 <- run_models(c(reg.formula1, reg.formula2)) 
mod.cl1 <- run_coef(mod1)

##DV means
means_dv<-c(rep(mean(dat$greens_voters1980, na.rm=TRUE),2), rep(mean(dat$greens_voters1983, na.rm=TRUE),2), 
            rep(mean(dat$greens_voters1987, na.rm=TRUE),2), rep(mean(dat$greens_voters1990, na.rm=TRUE),2))
sds_dv<-c(rep(sd(dat$greens_voters1980, na.rm=TRUE),2), rep(sd(dat$greens_voters1983, na.rm=TRUE),2), 
          rep(sd(dat$greens_voters1987, na.rm=TRUE),2), rep(sd(dat$greens_voters1990, na.rm=TRUE),2))


tableGreens <- texreg(
  l = list(mod1[[1]], mod1[[5]], mod1[[2]],  mod1[[6]], mod1[[3]], mod1[[7]], mod1[[4]], mod1[[8]]),
  override.se=list(mod.cl1[[1]][,2], mod.cl1[[5]][,2],  mod.cl1[[2]][,2], mod.cl1[[6]][,2],  mod.cl1[[3]][,2], mod.cl1[[7]][,2], mod.cl1[[4]][,2], mod.cl1[[8]][,2]),
  override.pval=list(mod.cl1[[1]][,4], mod.cl1[[5]][,4],  mod.cl1[[2]][,4], mod.cl1[[6]][,4], mod.cl1[[3]][,4],mod.cl1[[7]][,4],  mod.cl1[[4]][,4],mod.cl1[[8]][,4]),     
  custom.model.names=c("1980", "1980", "1983","1983", "1987","1987", "1990", "1990"), 
  custom.coef.names=c("Delta diversity", "Religious fractionalization (1939)", "Share expellees","Share in agriculture (1939)", "Share female (1939)",  
                      "Ln distance to Eastern border", "Distance to station (1950)", "Population density (1939)","Destruction", "Ln population (1939)"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:greens",
  omit.coef = "kreis_kenn1950_alt|Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE, 
  custom.gof.rows = list("DV Mean" = means_dv, "DV SD" = sds_dv)
)

##############################################
### Right panel of Figure 5 (Greens) #########
##############################################

plot.outGreens <- sapply(list(mod1[[5]], mod1[[6]], mod1[[7]], mod1[[8]]), function(x){
  coeftest(x, vcov.=hc1)[2,1:2] 
})

plot.outGreens <- as.data.frame(t(plot.outGreens))
plot.outGreens$DV<-c("1980", "1983", "1987", "1990")
plot.outGreens$IV<-rep("Delta diversity", 4)
plot.outGreens$lb <- plot.outGreens$Estimate - 1.96*plot.outGreens$`Std. Error`
plot.outGreens$ub <- plot.outGreens$Estimate + 1.96*plot.outGreens$`Std. Error`


figureGreens<-ggplot(plot.outGreens, aes(x = DV, y = Estimate, color = IV)) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_pointrange(aes(ymin = lb, ymax = ub), 
                  alpha = 0.75, 
                  position = position_dodge(width = 0.1), 
                  size = 0.5) +  # Set size to 0.5 for smaller circles
  scale_y_continuous(name = "Average marginal effect: Green Party vote") +
  # scale_x_discrete(name = "") +
  scale_color_manual(values = "black") +  # Set symbol color to black
  labs(color = "IV") +
  theme_bw() +
  theme(
    text=element_text(size=12, family="sans", color="black"),
    legend.position = "bottom",
    legend.text = element_text(size = 12,family="sans"),
    axis.text.x = element_text(angle = 47, hjust = 1, size=12, family="sans", color="black"),
    axis.title.x = element_blank(),  # Remove x-axis label
    legend.title = element_blank(),
    plot.margin = unit(c(.5,.5,.5,1), "cm")
  )


ggsave("images/Figure5right.tiff",figureGreens, width = 8, height = 6, dpi = 300, device = "tiff")

#####################################################################################
############# Table A18, panel A: Greens without county FEs OR with diocese FEs #####
#####################################################################################
dvGreens<-list("lm(greens_voters1980~", "lm(greens_voters1983~", "lm(greens_voters1987~", "lm(greens_voters1990~")

expvars1 <- list("DeltaFrac+ relfrac1939+srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939, data=dat)")

expvars2 <- list("DeltaFrac+ relfrac1939+srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(diocese), data=dat)")

reg.formula1 <- formulas(dvGreens, expvars1)
reg.formula2 <- formulas(dvGreens, expvars2)

mod1 <- run_models(c(reg.formula1, reg.formula2)) 
mod.cl1 <- run_coef(mod1)


tableGreensAlt1<- texreg(
  l = list(mod1[[1]], mod1[[5]], mod1[[2]],  mod1[[6]], mod1[[3]], mod1[[7]], mod1[[4]], mod1[[8]]),
  override.se=list(mod.cl1[[1]][,2], mod.cl1[[5]][,2],  mod.cl1[[2]][,2], mod.cl1[[6]][,2],  mod.cl1[[3]][,2], mod.cl1[[7]][,2], mod.cl1[[4]][,2], mod.cl1[[8]][,2]),
  override.pval=list(mod.cl1[[1]][,4], mod.cl1[[5]][,4],  mod.cl1[[2]][,4], mod.cl1[[6]][,4], mod.cl1[[3]][,4],mod.cl1[[7]][,4],  mod.cl1[[4]][,4],mod.cl1[[8]][,4]),     
  custom.model.names=c("1980", "1980", "1983","1983", "1987","1987", "1990", "1990"), 
  custom.coef.names=c("Delta diversity", "Religious fractionalization (1939)", "Share expellees", "Share in agriculture (1939)", "Share female (1939)",  
                      "Ln distance to Eastern border", "Distance to station (1950)", "Population density (1939)","Destruction", "Ln population (1939)"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:GreensAlt1",
  omit.coef = "kreis_kenn1950_alt|diocese|Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE 
)

######################################################################
############# Table A18, panel B: Greens without expellee share  #####
######################################################################

dvGreens<-list("lm(greens_voters1980~", "lm(greens_voters1983~", "lm(greens_voters1987~", "lm(greens_voters1990~")

expvars1 <- list("DeltaFrac+ relfrac1939+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939, data=dat)")

expvars2 <- list("DeltaFrac+ relfrac1939+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)")

reg.formula1 <- formulas(dvGreens, expvars1)
reg.formula2 <- formulas(dvGreens, expvars2)

mod1 <- run_models(c(reg.formula1, reg.formula2)) 
mod.cl1 <- run_coef(mod1)


tableGreensAlt2<- texreg(
  l = list(mod1[[1]], mod1[[5]], mod1[[2]],  mod1[[6]], mod1[[3]], mod1[[7]], mod1[[4]], mod1[[8]]),
  override.se=list(mod.cl1[[1]][,2], mod.cl1[[5]][,2],  mod.cl1[[2]][,2], mod.cl1[[6]][,2],  mod.cl1[[3]][,2], mod.cl1[[7]][,2], mod.cl1[[4]][,2], mod.cl1[[8]][,2]),
  override.pval=list(mod.cl1[[1]][,4], mod.cl1[[5]][,4],  mod.cl1[[2]][,4], mod.cl1[[6]][,4], mod.cl1[[3]][,4],mod.cl1[[7]][,4],  mod.cl1[[4]][,4],mod.cl1[[8]][,4]),     
  custom.model.names=c("1980", "1980", "1983","1983", "1987","1987", "1990", "1990"), 
  custom.coef.names=c("Delta diversity", "Religious fractionalization (1939)", "Share in agriculture (1939)", "Share female (1939)", "Ln Distance to Eastern Border", "Distance to station (1950)", "Population density (1939)","Destruction", "Ln population (1939)"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:selfemp",
  omit.coef = "kreis_kenn1950_alt|diocese|Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE, 
)

############################################################################
######## APPENDIX TABLE A19: HETEROGENEOUS EFFECTS by PREWAR DIVERSITY   ###
############################################################################

#New variables to average over outcomes for ease of presentation
dat$greens_voters1980s<-(dat$greens_voters1980+dat$greens_voters1983+dat$greens_voters1987)/3
dat$csu_voters_avg1980s<-(dat$csu_voters1980+dat$csu_voters1983+dat$csu_voters1987)/3

dvsHE<-list("lm(attendants_rel1970s~", "lm(preligious_70s~", "lm(pdivorce_tomar_70~", "lm(greens_voters1980s~", "lm(csu_voters_avg1980s~")
expvarsHE2 <- list("DeltaFrac+relfrac1939*DeltaFrac+relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)")

reg.formulaHE <- formulas(dvsHE, expvarsHE2)

modHE <- run_models(reg.formulaHE) 
mod.clHE <- run_coef(modHE)

tableHE2 <- texreg(
  l = list(modHE[[1]], modHE[[2]], modHE[[3]],  modHE[[4]], modHE[[5]]),
  override.se=list(mod.clHE[[1]][,2], mod.clHE[[2]][,2],  mod.clHE[[3]][,2], mod.clHE[[4]][,2],  mod.clHE[[5]][,2]),
  override.pval=list(mod.clHE[[1]][,4], mod.clHE[[2]][,4],  mod.clHE[[3]][,4], mod.clHE[[4]][,4], mod.clHE[[5]][,4]),     
  custom.model.names = c("mass attendantce 1970", "church affiliation 1970", "Divorces 1970", "Greens Voteshare 1980s", "CSU voteshare 1980s"),
  custom.coef.names=c("Delta diversity",  "Rel frac (1939)", "Share refugees",  "Share in agriculture (1939)","Share female (1939)",  
                      "Ln distance to Eastern border", "Distance to station (1950)", "Population density (1939)","Destruction" ,"Ln population (1939)", "Rel frac (1939) * Delta diversity"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:heterogeneity.diversity",
  omit.coef = "kreis_kenn1950_alt|Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE
)


#######################################################################################################
###### Figure A3: Plots of interaction effects using Interflex package  (Hainmueller, Mummolo, Xu) ####
#######################################################################################################


theme_set(theme_bw()) # for white background plots

moderator <- "relfrac1939"
controlsHE <- c("srefugees1950", "sagri1939", "sfemale1939", "lndistBorderEast", "diststation1950","popdens1939", "dest1945", "lnpop1939")

set.seed(1234) #because we are using nonparametric sampling with "kernel" we need to ensure that results remain the same across runs.

#jpeg("images/inter_attend_frac2.jpeg", width = 800, height = 600)  # Set file name and dimensions
tiff("images/FigureA3a.tiff", width = 800, height = 600)  
interflex(estimator = "kernel", 
          data = dat,
          Y = "attendants_rel1970s", 
          D = "DeltaFrac", 
          X = moderator, 
          treat.type = "continuous", 
          Z = controlsHE, 
          FE = "kreis_kenn1950_alt", 
          full.moderate = FALSE,
          weights = NULL, 
          na.rm = TRUE, 
          CI = TRUE, 
          Ylabel = "Religious Affiliation 1970", 
          Dlabel = "Delta diversity", 
          Xlabel="Religious fractionalization in 1939", 
          main = "Outcome: Mass attendance in 1970", 
          cex.main = 1.2,
)  
dev.off()


#jpeg("images/inter_affiliation_frac2.jpeg", width = 800, height = 600)  # Set file name and dimensions
tiff("images/FigureA3b.tiff", width = 800, height = 600)  
interflex(estimator = "kernel", 
          data = dat,
          Y = "preligious_70s", 
          D = "DeltaFrac", 
          X = moderator, 
          treat.type = "continuous", 
          Z = controlsHE, 
          FE = "kreis_kenn1950_alt", 
          full.moderate = FALSE,
          weights = NULL, 
          na.rm = TRUE, 
          CI = TRUE, 
          Ylabel = "Religious affiliation 1970", 
          Dlabel = "Delta diversity", 
          Xlabel="Religious fractionalization in 1939", 
          main = "Outcome: Religious affiliation in 1970", 
          cex.main = 1.2,
)  
dev.off()


#jpeg("images/inter_divorce_frac2.jpeg", width = 800, height = 600)  # Set file name and dimensions
tiff("images/FigureA3c.tiff", width = 800, height = 600)  
interflex(estimator = "kernel", 
          data = dat,
          Y = "pdivorce_tomar_70", 
          D = "DeltaFrac", 
          X = moderator, 
          treat.type = "continuous", 
          Z = controlsHE, 
          FE = "kreis_kenn1950_alt", 
          full.moderate = FALSE,
          weights = NULL, 
          na.rm = TRUE, 
          CI = TRUE, 
          Ylabel = "Divorce rates in 1970", 
          Dlabel = "Delta diversity", 
          Xlabel="Religious fractionalization in 1939", 
          main = "Outcome: Divorce rates in 1970", 
          cex.main = 1.2,
)    
dev.off()

#jpeg("images/inter_divorce87_frac2.jpeg", width = 800, height = 600)  # Set file name and dimensions
tiff("images/FigureA3d.tiff", width = 800, height = 600)  
interflex(estimator = "kernel", 
          data = dat,
          Y = "pdivorce_tomar_87", 
          D = "DeltaFrac", 
          X = moderator, 
          treat.type = "continuous", 
          Z = controlsHE, 
          FE = "kreis_kenn1950_alt", 
          full.moderate = FALSE,
          weights = NULL, 
          na.rm = TRUE, 
          CI = TRUE, 
          Ylabel = "Divorce rates in 1970", 
          Dlabel = "Delta diversity", 
          Xlabel="Religious fractionalization in 1939", 
          main = "Outcome: Divorce rates in 1987", 
          cex.main = 1.2,
)    
dev.off()


#jpeg("images/inter_Greens_frac2.jpeg", width = 800, height = 600)  # Set file name and dimensions
tiff("images/FigureA3e.tiff", width = 800, height = 600)  
interflex(estimator = "kernel", 
          data = dat,
          Y = "greens_voters1980s", 
          D = "DeltaFrac", 
          X = moderator, 
          treat.type = "continuous", 
          Z = controlsHE, 
          FE = "kreis_kenn1950_alt", 
          full.moderate = FALSE,
          weights = NULL, 
          na.rm = TRUE, 
          CI = TRUE, 
          Ylabel = "Green Party voteshare", 
          Dlabel = "Delta diversity", 
          Xlabel="Religious fractionalization in 1939", 
          main = "Outcome: Green Party voteshare (1980-1990 avg)", 
          cex.main = 1.2,
)    
dev.off()

#jpeg("images/inter_CSU_frac2.jpeg", width = 800, height = 600)  # Set file name and dimensions
tiff("images/FigureA3f.tiff", width = 800, height = 600)  
interflex(estimator = "kernel", 
          data = dat,
          Y = "csu_voters_avg1980s", 
          D = "DeltaFrac", 
          X = moderator, 
          treat.type = "continuous", 
          Z = controlsHE, 
          FE = "kreis_kenn1950_alt", 
          full.moderate = FALSE,
          weights = NULL, 
          na.rm = TRUE, 
          CI = TRUE, 
          Ylabel = "CSU voteshare", 
          Dlabel = "Delta diversity", 
          Xlabel="Religious fractionalization in 1939", 
          main = "Outcome: CSU voteshare (1980-1990 avg)", 
          cex.main = 1.2,
)
dev.off()
  
##########################################################################################
###### Appendix Section H, Table A20: SENSITIVITY ANALYSIS: Cinelli & Hazlett 2019 #######
##########################################################################################

#Note ovb_minimal_reporting() prints only one benchmark covariate, so to check results for both, you can comment out one of the benchmark_covariates options 

#1. For church attendance in 1970:

mod_churchatt<-lm(attendants_rel1970s~DeltaFrac+ srefugees1950+relfrac1939+ sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)

attendance.sensitivity <- sensemakr(model = mod_churchatt, 
                                    treatment = "DeltaFrac",
                                    benchmark_covariates = "sfemale1939", #to obtain bounds for sfemale1939, remove # from line 1393 and comment out line 1394 instead
                                    #benchmark_covariates = c("sagri1939", "sfemale1939"), #Figure below uses sensemakr output with two benchmark covariates, sagri1939 and sfemale1939
                                    kd = 1:3)
ovb_minimal_reporting(attendance.sensitivity, format = "latex")

jpeg("images/sensitivity_church_attendance2FE.jpeg", width = 600, height = 600)
plot(attendance.sensitivity, main="Church attendance (1970) sensitivity") 
dev.off()

#2. Affiliation with the church in 1970.
mod_affil<-lm(preligious_70s~DeltaFrac+ relfrac1939+srefugees1950+ sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)

affiliation.sensitivity <- sensemakr(model = mod_affil, 
                                     treatment = "DeltaFrac",
                                     benchmark_covariates = "sfemale1939", #to obtain bounds for sfemale1939, remove # from line 1407 and comment out line 1408 instead
                                     #benchmark_covariates = c("sagri1939","sfemale1939"),
                                     kd = 1:3)

ovb_minimal_reporting(affiliation.sensitivity, format = "latex")

#Figure uses sensemakr output with two benchmark covariates, sagri1939 and sfemale1939
jpeg("Laender/Vertriebene_Bayern/Analysis/Figures/sensitivity_affiliation2FE.jpeg", width = 600, height = 600)
plot(affiliation.sensitivity, main="Church affiliation (1970) sensitivity")
dev.off()

#3. Divorce rates in 1970 (divorces to marriages)

mod_div<-lm(pdivorce_tomar_70~DeltaFrac+ relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)

divorce.sensitivity <- sensemakr(model = mod_div, 
                                 treatment = "DeltaFrac",
                                 benchmark_covariates = "sfemale1939", #to obtain bounds for sfemale1939, remove # from line 1423 and comment out line 1424 instead
                                 #benchmark_covariates = c("sagri1939","sfemale1939"),
                                 kd = 1:3)

ovb_minimal_reporting(divorce.sensitivity, format = "latex")

#Figure uses sensemakr output with two benchmark covariates, sagri1939 and sfemale1939
jpeg("images/sensitivity_divorce2FE.jpeg", width = 600, height = 600)
plot(divorce.sensitivity, main="Divorce rates (1970) sensitivity")
dev.off()

#4-5. CSU & Greens vote (Averaging 1980, 1983, 1987)
dat$greens_voters1980s<-(dat$greens_voters1980+dat$greens_voters1983+dat$greens_voters1987)/3
dat$csu_voters_avg1980s<-(dat$csu_voters1980+dat$csu_voters1983+dat$csu_voters1987)/3

dvs<-list("lm(csu_voters_avg1980s~", "lm(greens_voters1980s~")
expvars <- list("DeltaFrac+relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)")
reg.formula.vote <- formulas(dvs, expvars)
mod.vote <- run_models(reg.formula.vote) 
mod.cl.vote <- run_coef(mod.vote)

csu.sensitivity <- sensemakr(model = mod.vote[[1]], 
                             treatment = "DeltaFrac",
                             benchmark_covariates = "sfemale1939", #to obtain bounds for sfemale1939, remove # from this line and comment out the next line
                             #benchmark_covariates = c("sagri1939","sfemale1939"),
                             kd = 1:3)

ovb_minimal_reporting(csu.sensitivity, format = "latex")

#Figure uses sensemakr output with two benchmark covariates, sagri1939 and sfemale1939
jpeg("images/sensitivity_CSU2FE.jpeg", width = 600, height = 600)
plot(csu.sensitivity, main="CSU vote (1980s) sensitivity")
dev.off()

greens.sensitivity <- sensemakr(model = mod.vote[[2]], 
                                treatment = "DeltaFrac",
                                benchmark_covariates = c("sfemale1939"), #to obtain bounds for sfemale1939, remove # from this line and comment out the next line
                                #benchmark_covariates = c("sagri1939","sfemale1939"),
                                kd = 1:3)

ovb_minimal_reporting(greens.sensitivity, format = "latex")

#Figure uses sensemakr output with two benchmark covariates, sagri1939 and sfemale1939
jpeg("images/sensitivity_Greens2FE.jpeg", width = 600, height = 600)
plot(greens.sensitivity, main="Green Party vote (1980s) sensitivity")
dev.off()

######################################################
##### Appendix Table A26: Economic indicators ########
######################################################


dvs<-list("lm(InAgric_70~", "lm(InManuf70~", "lm(HandelVerkehr70~", "lm(SelbstStandige_70~")

expvars2 <- list("DeltaFrac+ relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)")

reg.formulaAlt <- formulas(dvs, expvars2)

mod1 <- run_models(reg.formulaAlt) 
mod.cl1 <- run_coef(mod1)

means_dv<-c(mean(dat$InAgric_70, na.rm = TRUE), mean(dat$InManuf70, na.rm = TRUE), mean(dat$HandelVerkehr70, na.rm = TRUE), mean(dat$SelbstStandige_70, na.rm = TRUE))
sds_dv<-c(sd(dat$InAgric_70, na.rm = TRUE), sd(dat$InManuf70, na.rm = TRUE), sd(dat$HandelVerkehr70, na.rm = TRUE), sd(dat$SelbstStandige_70, na.rm = TRUE))


table_econalt <- texreg(
  l = list(mod1[[1]], mod1[[2]], mod1[[3]], mod1[[4]]),
  override.se=list(mod.cl1[[1]][,2], mod.cl1[[2]][,2], mod.cl1[[3]][,2],  mod.cl1[[4]][,2]),
  override.pval=list(mod.cl1[[1]][,4],  mod.cl1[[2]][,4], mod.cl1[[3]][,4], mod.cl1[[4]][,4]),     
  custom.model.names=c("In agriculture", "In manufacturing" , "In trade and services", "Self-employed"), 
  custom.coef.names=c("Delta diversity", "Religious fractionalization (1939)","Share expellees", "Share in agriculture (1939)", "Share female (1939)",  
                     "Ln distance to Eastern border", "Distance to station (1950)", "Population density (1939)", "Destruction", "Ln population (1939)"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  label = "tab:econ_dev",
  omit.coef = "kreis_kenn1950_alt|Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE, 
  custom.gof.rows = list("DV Mean" = means_dv, "DV SD" = sds_dv)
)


#################################################################################################
################### Table A27: ACDE when controlling for mediators (sequential g-estimation) ####
#################################################################################################

#Centering mediator variables: 
dat$centered_InAgric_70 <- dat$InAgric_70 - mean(dat$InAgric_70, na.rm = TRUE)
dat$centered_SelbstStandige_70 <- dat$SelbstStandige_70 - mean(dat$SelbstStandige_70, na.rm = TRUE)
dat$centered_HandelVerkehr70 <- dat$HandelVerkehr70 - mean(dat$HandelVerkehr70, na.rm = TRUE)
dat$centered_relfrac1970 <- dat$relfrac1970 - mean(dat$relfrac1970, na.rm = TRUE)

dat$centered_churchatt1970s <- dat$attendants_rel1970s - mean(dat$attendants_rel1970s, na.rm = TRUE)
dat$centered_churchatt1970SQ <- dat$centered_churchatt1970s^2

dat$centered_affiliation1970s <- dat$preligious_70s - mean(dat$preligious_70s, na.rm = TRUE)
dat$centered_affiliation1970SQ <- dat$centered_affiliation1970s^2


#Subset to municipalities in the Archdiocese of Munich-Freising, for which church attendance data exist
attend_subset<-subset(dat, !is.na(dat$centered_churchatt1970s))

#Values for table A26: DV: Mass attendance  in 1970

#Obtain baseline estimates (same as in the paper's main tables)
ate.mod <- lm(attendants_rel1970s ~ DeltaFrac+ relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)
round(summary(ate.mod)$coefficients["DeltaFrac", "Estimate"], 3) #Baseline estimate for mass attendance
round(summary(ate.mod)$coefficients["DeltaFrac", "Pr(>|t|)"], 3) #p-value indicated by the number of stars
round(confint(ate.mod)[2,1:2], 3) #95% CI of the baseline estimate

#Obtain results after fixing different mediators: yvar ~ treat + xvars | zvars | mvars
church_med_agr<-attendants_rel1970s~ DeltaFrac+ relfrac1939+srefugees1950+ sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt) |  pAuslaender1970  +SelbstStandige_70+HandelVerkehr70+relfrac1970 | centered_InAgric_70
church_med_handel<-attendants_rel1970s~ DeltaFrac+ relfrac1939+srefugees1950+ sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt) |  pAuslaender1970  +SelbstStandige_70+InAgric_70+relfrac1970 | centered_HandelVerkehr70
church_med_semp<-attendants_rel1970s~ DeltaFrac+relfrac1939+srefugees1950+  sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt) |  pAuslaender1970  +HandelVerkehr70+InAgric_70+relfrac1970 | centered_SelbstStandige_70

direct1 <- sequential_g(church_med_agr, data = dat)
direct2 <- sequential_g(church_med_handel, data = dat)
direct3 <- sequential_g(church_med_semp, data = dat)

# Collecting estimates, p-values, and confidence intervals
results <- data.frame(
  Model = c("Baseline", rep("Sequential g-estimate, ACDE", 3)),  
  Mediator = c("None", "Agriculture", "Trade", "Self-Employment"),
  Estimate = c(
    paste0(round(summary(ate.mod)$coefficients["DeltaFrac", "Estimate"], 3), pval_stars(summary(ate.mod)$coefficients["DeltaFrac", "Pr(>|t|)"])),
    paste0(round(summary(direct1)$coefficients["DeltaFrac", "Estimate"], 3), pval_stars(summary(direct1)$coefficients["DeltaFrac", "Pr(>|t|)"])),
    paste0(round(summary(direct2)$coefficients["DeltaFrac", "Estimate"], 3), pval_stars(summary(direct2)$coefficients["DeltaFrac", "Pr(>|t|)"])),
    paste0(round(summary(direct3)$coefficients["DeltaFrac", "Estimate"], 3), pval_stars(summary(direct3)$coefficients["DeltaFrac", "Pr(>|t|)"]))
  ),
  `95% CI` = c(
    paste0("[", round(confint(ate.mod)[2, 1], 3), ", ", round(confint(ate.mod)[2, 2], 3), "]"),
    paste0("[", round(confint(direct1)[2, 1], 3), ", ", round(confint(direct1)[2, 2], 3), "]"),
    paste0("[", round(confint(direct2)[2, 1], 3), ", ", round(confint(direct2)[2, 2], 3), "]"),
    paste0("[", round(confint(direct3)[2, 1], 3), ", ", round(confint(direct3)[2, 2], 3), "]")  
  )
)

xtable(results, caption = "DV: Mass attendance in 1970")

## Row 2: DV: Divorce rates in 1970.

#baseline estimation:
ate.mod2 <- lm(pdivorce_tomar_70 ~ DeltaFrac+ relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)

#Obtain results after fixing different mediators: yvar ~ treat + xvars | zvars | mvars

divorce_med_agr<-pdivorce_tomar_70~ DeltaFrac+ relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt) |  pAuslaender1970  +SelbstStandige_70+HandelVerkehr70+relfrac1970 | centered_InAgric_70
divorce_med_handel<-pdivorce_tomar_70~ DeltaFrac+ relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt) |  pAuslaender1970  +SelbstStandige_70+InAgric_70+relfrac1970 | centered_HandelVerkehr70
divorce_med_semp<-pdivorce_tomar_70~ DeltaFrac+ relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt) |  pAuslaender1970  +HandelVerkehr70+InAgric_70+relfrac1970 | centered_SelbstStandige_70

direct4 <- sequential_g(divorce_med_agr, data = dat)

direct5 <- sequential_g(divorce_med_handel, data = dat)

direct6 <- sequential_g(divorce_med_semp, data = dat)

results2 <- data.frame(
  Model = c("Baseline", rep("Sequential g-estimate, ACDE", 3)),  
  Mediator = c("None", "Agriculture", "Trade", "Self-Employment"),
  Estimate = c(
    paste0(round(summary(ate.mod2)$coefficients["DeltaFrac", "Estimate"], 3), pval_stars(summary(ate.mod2)$coefficients["DeltaFrac", "Pr(>|t|)"])),
    paste0(round(summary(direct4)$coefficients["DeltaFrac", "Estimate"], 3), pval_stars(summary(direct4)$coefficients["DeltaFrac", "Pr(>|t|)"])),
    paste0(round(summary(direct5)$coefficients["DeltaFrac", "Estimate"], 3), pval_stars(summary(direct5)$coefficients["DeltaFrac", "Pr(>|t|)"])),
    paste0(round(summary(direct6)$coefficients["DeltaFrac", "Estimate"], 3), pval_stars(summary(direct6)$coefficients["DeltaFrac", "Pr(>|t|)"]))
  ),
  `95% CI` = c(
    paste0("[", round(confint(ate.mod2)[2, 1], 3), ", ", round(confint(ate.mod2)[2, 2], 3), "]"),
    paste0("[", round(confint(direct4)[2, 1], 3), ", ", round(confint(direct4)[2, 2], 3), "]"),
    paste0("[", round(confint(direct5)[2, 1], 3), ", ", round(confint(direct5)[2, 2], 3), "]"),
    paste0("[", round(confint(direct6)[2, 1], 3), ", ", round(confint(direct6)[2, 2], 3), "]")  
  )
)

xtable(results2, caption = "DV: Divorce rates in 1970")

## Row 3: DV: Church affiliation in 1970: 
#baseline estimate
ate.mod3 <- lm(preligious_70s ~ DeltaFrac+ relfrac1939+ srefugees1950+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)

#Obtain results after fixing different mediators: yvar ~ treat + xvars | zvars | mvars

affil_med_agr<-preligious_70s~ DeltaFrac+ relfrac1939+srefugees1950+ sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt) |  pAuslaender1970  +SelbstStandige_70+HandelVerkehr70+relfrac1970 | centered_InAgric_70
affil_med_handel<-preligious_70s~ DeltaFrac+ relfrac1939+srefugees1950+ sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt) |  pAuslaender1970  +SelbstStandige_70+InAgric_70+relfrac1970 | centered_HandelVerkehr70
affil_med_semp<-preligious_70s~ DeltaFrac+ relfrac1939+srefugees1950+ sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt) |  pAuslaender1970  +HandelVerkehr70+InAgric_70+relfrac1970 | centered_SelbstStandige_70


direct7 <- sequential_g(affil_med_agr, data = dat)
direct8 <- sequential_g(affil_med_handel, data = dat)
direct9 <- sequential_g(affil_med_semp, data = dat)

results3 <- data.frame(
  Model = c("Baseline", rep("Sequential g-estimate, ACDE", 3)),  
  Mediator = c("None", "Agriculture", "Trade", "Self-Employment"),
  Estimate = c(
    paste0(round(summary(ate.mod3)$coefficients["DeltaFrac", "Estimate"], 3), pval_stars(summary(ate.mod3)$coefficients["DeltaFrac", "Pr(>|t|)"])),
    paste0(round(summary(direct7)$coefficients["DeltaFrac", "Estimate"], 3), pval_stars(summary(direct7)$coefficients["DeltaFrac", "Pr(>|t|)"])),
    paste0(round(summary(direct8)$coefficients["DeltaFrac", "Estimate"], 3), pval_stars(summary(direct8)$coefficients["DeltaFrac", "Pr(>|t|)"])),
    paste0(round(summary(direct9)$coefficients["DeltaFrac", "Estimate"], 3), pval_stars(summary(direct9)$coefficients["DeltaFrac", "Pr(>|t|)"]))
  ),
  `95% CI` = c(
    paste0("[", round(confint(ate.mod3)[2, 1], 3), ", ", round(confint(ate.mod3)[2, 2], 3), "]"),
    paste0("[", round(confint(direct7)[2, 1], 3), ", ", round(confint(direct7)[2, 2], 3), "]"),
    paste0("[", round(confint(direct8)[2, 1], 3), ", ", round(confint(direct8)[2, 2], 3), "]"),
    paste0("[", round(confint(direct9)[2, 1], 3), ", ", round(confint(direct9)[2, 2], 3), "]")  
  )
)

xtable(results3, caption = "DV: Church Affiliation in 1970")

#Figure A.5: Sensitivity for share in agriculture as mediator. 
direct1 <- sequential_g(attendants_rel1970s~ DeltaFrac+srefugees1950+ relfrac1939+ sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1950 |  pAuslaender1970  +SelbstStandige_70+HandelVerkehr70+relfrac1970 | centered_InAgric_70,
                        dat)
out_sens <- cdesens(direct1, var = "DeltaFrac")

plot(out_sens)

###############################################################################################
############ Table A28: ACDE of Share expellees   ####
###############################################################################################
#Note these models use lnpop1939 (pre-treatment measure)

dat$csu_voters_avg1950s<-(dat$csu_voters1953+dat$csu_voters1957)/3

dat$DeltaFracSQ<-dat$DeltaFrac^2 #squared term for delta fraq

#Step 1: Obtain baseline estimates for treatment Share expellees (srefugees1950) [no DeltaFrac here]

ate.mod_mass <- lm(attendants_rel1970s ~ srefugees1950+relfrac1939+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)
#The effect of refugee presence is negative and significant.
ate.mod_divorce <- lm(pdivorce_tomar_70 ~ srefugees1950+relfrac1939+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)

ate.mod_affiliation <- lm(preligious_70s ~ srefugees1950+relfrac1939+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)
ate.mod_csu <- lm(csu_voters_avg1950s ~ srefugees1950+relfrac1939+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)

ate.mod_greens <- lm(greens_voters1980s ~ srefugees1950+relfrac1939+sagri1939+sfemale1939 +lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt), data=dat)


#Step 2: estimate controlled direct effect of refugees fixing the value of the medicator, centered_DeltaFrac, at the mean

#Structure: yvar ~ treat + xvars | zvars | mvars
#Our intermediate covariate for 1950 are per capita municipal revenue, share of the population in agriculture, share of females, share of the elderly (abover 65), and share commuters:
formula_mass<-attendants_rel1970s~ srefugees1950+relfrac1939+ sagri1939 +sfemale1939+lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt) |  
  taxRevenuePC1950+InAgric_50+sfemale1950+PT65andolder1950+ShareCommuters1950 | DeltaFrac + DeltaFracSQ

formula_divorce<-pdivorce_tomar_70~ srefugees1950+relfrac1939+ sagri1939 +sfemale1939+lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt) |  
  taxRevenuePC1950+InAgric_50+sfemale1950+PT65andolder1950+ShareCommuters1950 | DeltaFrac + DeltaFracSQ

formula_affiliation<-preligious_70s~ srefugees1950+relfrac1939+ sagri1939 +sfemale1939+lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt) |  
  taxRevenuePC1950+InAgric_50+sfemale1950+PT65andolder1950+ShareCommuters1950 | DeltaFrac + DeltaFracSQ

formula_csu<-csu_voters_avg1950s~ srefugees1950+relfrac1939+ sagri1939 +sfemale1939+lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt) |  
  taxRevenuePC1950+InAgric_50+sfemale1950+PT65andolder1950+ShareCommuters1950 | DeltaFrac + DeltaFracSQ

formula_greens<-greens_voters1980s~ srefugees1950+relfrac1939+ sagri1939 +sfemale1939+lndistBorderEast+diststation1950+popdens1939+dest1945+lnpop1939+as.factor(kreis_kenn1950_alt) |  
  taxRevenuePC1950+InAgric_50+sfemale1950+PT65andolder1950+ShareCommuters1950 | DeltaFrac + DeltaFracSQ


direct1 <- sequential_g(formula_mass, data = dat)
direct2 <- sequential_g(formula_divorce, data = dat)
direct3 <- sequential_g(formula_affiliation, data = dat)
direct4 <- sequential_g(formula_csu, data = dat)
direct5 <- sequential_g(formula_greens, data = dat)

# Collecting estimates, p-values, and confidence intervals
resultsA28 <- data.frame(
  Model = rep(c("Baseline","Sequential g-estimate, ACDE"), 5), 
  Mediator = c("None", "Delta Diversity", "None", "Delta Diversity", "None", "Delta Diversity","None", "Delta Diversity", "None", "Delta Diversity"),
  Estimate = c(
    paste0(round(summary(ate.mod_mass)$coefficients["srefugees1950", "Estimate"], 3), pval_stars(summary(ate.mod_mass)$coefficients["srefugees1950", "Pr(>|t|)"])),
    paste0(round(summary(direct1)$coefficients["srefugees1950", "Estimate"], 3), pval_stars(summary(direct1)$coefficients["srefugees1950", "Pr(>|t|)"])),
    paste0(round(summary(ate.mod_divorce)$coefficients["srefugees1950", "Estimate"], 3), pval_stars(summary(ate.mod_divorce)$coefficients["srefugees1950", "Pr(>|t|)"])),
    paste0(round(summary(direct2)$coefficients["srefugees1950", "Estimate"], 3), pval_stars(summary(direct2)$coefficients["srefugees1950", "Pr(>|t|)"])),
    paste0(round(summary(ate.mod_affiliation)$coefficients["srefugees1950", "Estimate"], 3), pval_stars(summary(ate.mod_affiliation)$coefficients["srefugees1950", "Pr(>|t|)"])),
    paste0(round(summary(direct3)$coefficients["srefugees1950", "Estimate"], 3), pval_stars(summary(direct3)$coefficients["srefugees1950", "Pr(>|t|)"])),
    paste0(round(summary(ate.mod_csu)$coefficients["srefugees1950", "Estimate"], 3), pval_stars(summary(ate.mod_csu)$coefficients["srefugees1950", "Pr(>|t|)"])),
    paste0(round(summary(direct4)$coefficients["srefugees1950", "Estimate"], 3), pval_stars(summary(direct4)$coefficients["srefugees1950", "Pr(>|t|)"])),
    paste0(round(summary(ate.mod_greens)$coefficients["srefugees1950", "Estimate"], 3), pval_stars(summary(ate.mod_greens)$coefficients["srefugees1950", "Pr(>|t|)"])),
    paste0(round(summary(direct5)$coefficients["srefugees1950", "Estimate"], 3), pval_stars(summary(direct5)$coefficients["srefugees1950", "Pr(>|t|)"]))
  ),
  `95% CI` = c(
    paste0("[", round(confint(ate.mod_mass)[2, 1], 3), ", ", round(confint(ate.mod_mass)[2, 2], 3), "]"),
    paste0("[", round(confint(direct1)[2, 1], 3), ", ", round(confint(direct1)[2, 2], 3), "]"),
    paste0("[", round(confint(ate.mod_divorce)[2, 1], 3), ", ", round(confint(ate.mod_divorce)[2, 2], 3), "]"),
    paste0("[", round(confint(direct2)[2, 1], 3), ", ", round(confint(direct2)[2, 2], 3), "]"),
    paste0("[", round(confint(ate.mod_affiliation)[2, 1], 3), ", ", round(confint(ate.mod_affiliation)[2, 2], 3), "]"),
    paste0("[", round(confint(direct3)[2, 1], 3), ", ", round(confint(direct3)[2, 2], 3), "]")  ,
    paste0("[", round(confint(ate.mod_csu)[2, 1], 3), ", ", round(confint(ate.mod_csu)[2, 2], 3), "]"),
    paste0("[", round(confint(direct4)[2, 1], 3), ", ", round(confint(direct4)[2, 2], 3), "]"),
    paste0("[", round(confint(ate.mod_greens)[2, 1], 3), ", ", round(confint(ate.mod_greens)[2, 2], 3), "]"),
    paste0("[", round(confint(direct5)[2, 1], 3), ", ", round(confint(direct5)[2, 2], 3), "]")
    
  )
)

xtable(resultsA28, caption = "ACDE of Share expellees when controlling for post-treatment confounding")


############################################################
######## Table A29, Panel A: County-level analyses (R2) ####
############################################################

#create averages of electoral outcomes
dat$csu_voters_avg1950s<-(dat$csu_voters1953+dat$csu_voters1957)/3
dat$greens_voters1980s<-(dat$greens_voters1980+dat$greens_voters1983+dat$greens_voters1987)/3

#collapse data to the level of prewar county (same as 1950 county)
collapsed_data <- dat %>%
  group_by(kreis_kenn1950_alt) %>%
  summarize(
    preligious_70s = mean(preligious_70s, na.rm = TRUE),
    attendants_rel1970s = mean(attendants_rel1970s, na.rm = TRUE),
    pdivorce_tomar_70 = mean(pdivorce_tomar_70, na.rm = TRUE),
    csu_voters1950s = mean(csu_voters_avg1950s, na.rm = TRUE),
    greens_voters1980s = mean(greens_voters1980s, na.rm = TRUE),
    srefugees1950 = mean(srefugees1950, na.rm = TRUE),
    DeltaFrac = mean(DeltaFrac, na.rm = TRUE),
    relfrac1939 = mean(relfrac1939, na.rm = TRUE),
    sagri1939 = mean(sagri1939, na.rm = TRUE),
    sfemale1939 = mean(sfemale1939, na.rm = TRUE),
    lndistBorderEast = log(mean(distBorderEast, na.rm = TRUE)),
    diststation1950 = mean(diststation1950, na.rm = TRUE),
    popdens1939 = mean(popdens1939, na.rm = TRUE), 
    dest1945=mean(dest1945, na.rm = TRUE), 
    lnpop1939=log(mean(pop1939, na.rm = TRUE)), 
    lnpop1950=log(mean(pop1950, na.rm = TRUE)), 
  )


dvs<-c("lm(attendants_rel1970s~", "lm(preligious_70s~", "lm(pdivorce_tomar_70~", "lm(csu_voters1950s~", "lm(greens_voters1980s~")
expvars <- list("DeltaFrac + relfrac1939  +  srefugees1950+sagri1939 + sfemale1939 + lndistBorderEast + diststation1950 + popdens1939+dest1945+lnpop1939, data=collapsed_data)")

reg.formula <- formulas(dvs, expvars)

county_mods <- run_models(reg.formula) 
county_mods.cl <- run_coef(county_mods)

tableA29panelA <- texreg(
  l = list(county_mods[[1]], county_mods[[2]], county_mods[[3]], county_mods[[4]],  county_mods[[5]]),
  override.se=list(county_mods.cl[[1]][,2], county_mods.cl[[2]][,2], county_mods.cl[[3]][,2],  county_mods.cl[[4]][,2], county_mods.cl[[5]][,2]),
  override.pval=list(county_mods.cl[[1]][,4],  county_mods.cl[[2]][,4], county_mods.cl[[3]][,4], county_mods.cl[[4]][,4], county_mods.cl[[5]][,4]),     
  custom.model.names = c("Church Attendance (1970)", "Religious Affiliation (1970)", "Divorce Rate (1970)", "CSU vote (1950s)", "Greens vote (1980s)"),
  custom.coef.names=c("Delta diversity", "Religious fractionalization (1939)",  "Share expellees", "Share in agriculture (1939)", "Share female (1939)",  
                      "Ln Distance to Eastern Border", "Distance to station (1950)", "Population density (1939)", "Destruction", "Ln population (1939)"),
  dcolumn = TRUE,
  booktabs = TRUE,
  stars = c(.01,.05,.1),
  digits=3,
  label = "tab:county_replication",
  omit.coef = "Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE
)

## Table A29, Panel B: Alternative specification without srefugees1950

dvs<-c("lm(attendants_rel1970s~", "lm(preligious_70s~", "lm(pdivorce_tomar_70~", "lm(csu_voters1950s~", "lm(greens_voters1980s~")
expvars <- list("DeltaFrac + relfrac1939  +  sagri1939 + sfemale1939 + lndistBorderEast + diststation1950 + popdens1939+dest1945+lnpop1939, data=collapsed_data)")

reg.formula <- formulas(dvs, expvars)

county_mods <- run_models(reg.formula) 
county_mods.cl <- run_coef(county_mods)

table.A29panelB <- texreg(
  l = list(county_mods[[1]], county_mods[[2]], county_mods[[3]], county_mods[[4]],  county_mods[[5]]),
  override.se=list(county_mods.cl[[1]][,2], county_mods.cl[[2]][,2], county_mods.cl[[3]][,2],  county_mods.cl[[4]][,2], county_mods.cl[[5]][,2]),
  override.pval=list(county_mods.cl[[1]][,4],  county_mods.cl[[2]][,4], county_mods.cl[[3]][,4], county_mods.cl[[4]][,4], county_mods.cl[[5]][,4]),     
  custom.model.names = c("Church Attendance (1970)", "Religious Affiliation (1970)", "Divorce Rate (1970)", "CSU vote (1950s)", "Greens vote (1980s)"),
  custom.coef.names=c("Delta diversity", "Religious fractionalization (1939)",   "Share in agriculture (1939)", "Share female (1939)",  
                      "Ln Distance to Eastern Border", "Distance to station (1950)", "Population density (1939)", "Destruction", "Ln population (1939)"),
  dcolumn = TRUE,
  booktabs = TRUE,
  digits=3,
  stars = c(.01,.05,.1),
  label = "tab:county_replication",
  omit.coef = "Intercept",
  include.rsquared = FALSE, 
  include.rmse=FALSE
)


